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CHAPTER  I 


PREFACE  TO  FIRST  VERSION,  VERSION  3.10 


Although  this  Is  the  first  QCPE  I3M-PC  copy  of  MOPAC  released,  it 
has  been  given  the  version-number  3.1.  This  Is  to  allow  users  to 
compare  results  with  other  copies  of  MOPAC  running  on  other  computers. 
The  results  of  all  version  3.1  copies  of  MOPAC,  be  they  from  IBM 
mainframes,  CRAY-XMP  supercomputers,  >r  VAX  1 1 -780 1 s  should  all  be 
virtually  Identical.  Even  though  IBM-PC  microcomputers  are  relatively 
Inexpensive,  the  results  of  MOPAC  calculations  are  as  valid  as  if  a 
CRAY-XMP  48  had  been  used.  In  other  words,  the  fact  that  an  IBM-PC  Is 
the  host  computer  does  not  mean  that  serious  research  cannot  be 
undertaken  using  MOPAC.  Similarly,  MNDO,  AMI  and  MIND0/3  results, 
published  In  professional  journals,  can  be  readily  verified,  provided 
the  systems  Involved  are  small  enough  to  be  run  on  the  PC  (limit:  7 
heavy  atoms  (non-hydrogens)  and  7  hydrogen  atoms). 

Caution  to  new  MOPAC  users 


Output  from  MOPAC  Is  only  as  good  as  the  input.  For  example.  If  a 
calculated  heat  of  formation  Is  significantly  different  from  that 
published,  most  likely  an  error  has  been  made  In  the  data  supplied  to 
MOPAC.  The  program  Is  research-quality,  and  significant  work  can  be 
done  using  It.  -  However,  painstaking  care  is  required  before 
publication-quality  results  can  be  obtained. 
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CHAPTER  1 
INTRODUCTION 


MOPAC  Is  a  general-purpose  serai-empirical  molecular  orbital  package 
for  the  study  of  chemical  reactions.  The  semi-empirical  Hamiltonians 
MNDO,  MINDO/3,  and  AMI  are  Implemented,  and  calculations  of  vibrational 
spectra,  thermodynamic  quantities,  isotopic  substitution  effects  and 
force  constants  for  molecules.  Ions,  and  radicals  are  combined  In  a 
fully  Integrated  package.  One  limitation  Is  that  RHF  radicals  cannot  be 
used  with  FORCE  In  the  default  set  of  options.  However,  users  can  run 
RHF  radicals  with  FORCE  if  the  maximum  number  of  allowed  atoms  Is 
reduced.  Within  the  electronic  part  of  the  calculation,  eigenvectors 
and  localized  orbitals,  chemical  bond  Indices,  charges,  molecular 
orbitals  etc.  are  calculable.  For  studying  chemical  reactions,  a 
transition-state  location  routine  and  a  transition  state  optimizing 
routine  are  available.  For  potential  users  to  get  the  most  out  of  the 
program  they  must  understand  how  the  program  works,  how  to  enter  data, 
how  to  Interpret  the  results,  and  what  to  do  when  things  go  wrong. 

While  MOPAC  calls  upon  many  concepts  In  quantum  theory  and 
thermodynamics  and  uses  some  fairly  advanced  mathematics,  the  user  who 
Is  not  familiar  with  these  specialized  topics  should  not  feel  excluded 
from  using  It.  On  the  contrary,  MOPAC  Is  written  with  the 
non-theoretician  In  mind.  To  this  end  the  data  are  kept  as  simple  as 
possible;  this  means  that  users  can  give  their  attention  to  the 
chemistry  Involved,  and  not  concern  themselves  with  quantum  and 
thermodynamic  exotica. 
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INTRODUCTION 


1.1  SUMMARY  OF  MOPAC  CAPABILITIES 

1.  MNDO,  MINDO/3,  and  AMI  Hamiltonians. 

2.  RHF  and  UHF  methods. 

3.  Limited  Configuration  Interaction 

1 .  4  configurations 

2.  Singlets,  Doublets,  and  Triplets 

3.  Excited  states 

4.  Geometry  optimizations,  etc.,  on  specified  states 

4.  Single  SCF  calculation 

5.  Geometry  optimization 

6.  Gradient  minimization 

7.  Transition  state  location 

8.  Reaction  path  coordinate  calculation 

9.  Force  constant  calculation 

10.  Normal  coordinate  analysis 

11.  Transition  dipole  calculation 

12.  Thermodynamic  properties  calculation 

13.  Localized  orbitals 

14.  Covalent  bond  orders 
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1.2  LAYOUT  OF  DATA 

This  section  Is  aimed  at  the  complete  novice  —  someone  who  knows 
nothing  at  all  about  the  structure  of  a  MOPAC  data-flle. 

First  of  all,  there  are  at  most  four  possible  types  of  data-flles 
for  MOPAC,  but  the  simplest  data-flle  Is  the  most  commonly  used.  Rather 
than  define  It,  an  example  Is  shown  below. 


1.2.1  Example  Of  Data  For  Ethylene 


Line 

1 

UHF  PULAY  MIND03 

VECTORS  DENSITY  LOCAL  T-300 

Line 

2 

EXAMPLE  OF  DATA 

FOR  MOPAC 

Line 

3 

MINDO/3  UHF  CLOSED-SHELL 

D2D 

ETHYLENE 

Line 

4a 

C 

Line 

4b 

C 

1.400118  1 

Line 

4C 

H 

1.098326  1 

123.572063 

1 

Line 

4d 

H 

1.098326  1 

123.572063 

1 

180.000000 

0 

2 

1 

3 

Line 

4e 

H 

1.098326  1 

123.572063 

1 

90.000000 

0 

1 

2 

3 

Line 

4f 

H 

1.098326  1 

123.572063 

1 

270.000000 

0 

1 

2 

3 

Line 

5 

As  can  be  seen,  the  first  three  lines  are  textual.  The  first  line 
consists  of  key-words  (here  seven  key-words  are  shown).  These  control 
the  calculation.  The  next  two  lines  are  comments  or  titles.  The  user 
might  want  to  put  the  name  of  the  molecule  and  why  It  Is  being  run  on 
these  two  lines. 

These  three  lines  are  obligatory.  If  no  name  or  comment  Is  wanted, 
leave  blank  lines.  If  no  key-words  are  specified,  leave  a  blank  line. 
A  common  error  Is  to  have  a  blank  line  before  the  key-word  line:  this 
error  Is  quite  tricky  to  find,  so  be  careful  not  to  have  four  lines 
before  the  start  of  the  geometric  data.  Whatever  Is  decided,  the  three 
lines,  blank  or  otherwise,  are  obligatory. 

The  next  set  of  lines  defines  the  geometry.  In  the  example,  the 
numbers  are  all  neatly  lined  up;  this  Is  not  necessary,  but  does  make  It 
easier  when  looking  for  errors  In  the  data.  The  geometry  Is  defined  In 
lines  4a  to  4f;  line  5  terminates  both  the  geometry  and  the  data-flle. 
Any  additional  data,  for  example  symmetry  data,  would  follow  line  5.  An 
explanation  of  the  geometry  definition  Is  given  In  the  chapter  "GEOMETRY 
SPECIFICATION". 
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Summarizing,  then,  the  structure  for  a  MOPAC  data-flle  Is: 


Line  1:  Key-Words.  (See  chapter  2  on  definitions  of  key-words) 

Line  2:  Title  of  the  calculation,  e.g.  the  name  of  the  molecule  or  Ion. 

Line  3:  Other  Information  describing  the  calculation. 

Line  4:  Internal  or  cartesian  coordinates  (See  chapter  on 
specification  of  geometry) 

Line  5:  Blank  line  to  terminate  the  geometry  definition. 

Other  layouts  for  data-f lies  Involve  additions  to  the  simple 
layout.  These  additions  occur  at  the  end  of  the  data-flle,  after  line 
5.  The  three  most  common  additions  are: 

(a)  Symmetry  data:  This  follows  the  geometric  data,  and  Is  enaed 
by  a  blank  line. 

(b)  Reaction  path:  After  all  geometry  and  symmetry  data  (If  any) 
are  read  in,  points  on  the  reaction  coordinate  are  defined. 

(c)  Saddle  data:  A  complete  second  geometry  is  Input.  The  second 
geometry  follows  the  first  geometry  and  symmetry  data  (if  any) 


CHAPTER  2 


KEY-WORDS 


2.1  SPECIFICATION  OF  KEY-WORDS 

All  control  data  are  entered  In  the  form  of  key-words,  which  form 
the  first  line  of  a  data-flle.  The  order  In  which  key-words  appear  Is 
not  Important  although  they  must  be  separated  by  a  space.  Some 
key-words  can  be  abbreviated  (for  example  1  ELECTRON  can  be  entered  as 
1ELECT)  but  In  general  the  full  key-word  Is  preferred,  primarily  In 
order  to  more  clearly  document  the  calculation,  and  secondarily  to 
obviate  the  possibility  that  an  abbreviated  key-word  might  not  be 
recognized.  If  there  Is  Insufficient  space  In  the  first  line  for  all 
the  key-words  needed,  then  consider  abbreviating  the  longer  words.  One 
type  of  key-word,  those  ending  In  an  equal  sign,  such  as,  BAR-0.05,  may 
not  be  abbreviated,  and  the  full  word  needs  to  be  supplied. 

If  two  key-words  which  are  incompatible,  like  UHF  and  C.I.-2,  are 
supplied,  or  a  key-word  which  Is  Incompatible  with  the  species  supplied, 
for  Instance  TRIPLET  and  a  methyl  radical,  then  error  trapping  will 
normally  occur,  and  an  error  message  printed.  This  usually  takes  an 
Insignificant  time,  so  data  are  quickly  checked  for  obvious  errors. 
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2.2  FULL  LIST  OF  KEY-WORDS  USED  IN  MOPAC 


'<1 
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1  ELECTRON-  FINAL  ONE-ELECTRON  MATRIX  TO  BE  PRINTED 

OSCF  -  READ  IN  DATA,  THEN  STOP 

1SCF  -  DO  ONE  SCF  AND  THEN  STOP 

AMI  -  THE  AMI  HAMILTONIAN  TO  BE  USED 

BAR-n.n  -  REDUCE  BAR  LENGTH  BY  A  MAXIMUM  OF  n.n 

BIRADICAL-  SYSTEM  HAS  TWO  UNPAIRED  ELECTRONS 

BONDS  -  FINAL  BOND-ORDER  MATRIX  TO  BE  PRINTED 

C.I.-2  -  A  MULTI-ELECTRON  CONFIGURATION  INTERACTION  SPECIFIED 

CHARGE=n  -  CHARGE  ON  SYSTEM  *  n  (e.g.  NH4  =>  CHARGE- 1) 

COMPFG  -  PRINT  HEAT  OF  FORMATION  CALCULATED  IN  COMPFG 
CYCLES  -  PERFORM  MAXIMUM  NUMBER  OF  CYCLES  IN  NLLSQ 
DCART  -  PRINT  DETAILS  OF  WORKING  IN  DCART 
DEBUG  -  DEBUG  OPTION  TURNED  ON 
DEBUGPULAY  PRINT  DETAILS  OF  WORKING  IN  PULAY 
DENOUT  -  DENSITY  MATRIX  OUTPUT  (CHANNEL  10) 

DENSITY  -  FINAL  DENSITY  MATRIX  TO  BE  PRINTED 
DERIV  -  PRINT  PART  OF  WORKING  IN  DERIV 

DFORCE  -  FORCE  CALCULATION  SPECIFIED,  ALSO  PRINT  FORCE  MATRIX. 
DOUBLET  -  RHF  DOUBLET  STATE  REQUIRED 
DUMP-n  -  RESTART  FILES  TO  BE  WRITTEN  EVERY  n  SECONDS 
ECHO  -  DATA  ARE  ECHOED  BACK  BEFORE  CALCULATION  STARTS 

EIGS  -  ALL  EIGENVALUES  IN  ITER  TO  BE  PRINTED 

ENPART  -  ENERGY  TO  BE  PARTITIONED  INTO  COMPONENTS 
ESR  -  RHF  UNPAIREO  SPIN  DENSITY  TO  BE  CALCULATED 
EXCITED  -  FIRST  EXCITED  SINGLET  STATE  IS  TO  BE  OPTIMIZED 
FILL=n  -  IN  RHF  OPEN  AND  CLOSED  SHELL,  FORCE  M.O.  n 
TO  BE  FILLED 

FLEPO  -  PRINT  DETAILS  OF  GEOMETRY  OPTIMIZATION 

FMAT  -  PRINT  DETAILS  OF  WORKING  IN  FMAT 

FOCK  -  LAST  FOCK  MATRIX  TO  BE  PRINTED 

FORCE  -  FORCE  CALCULATION  SPECIFIED 

FULSCF  -  FULL  SCF  CALCN'S  TO  BE  DONE  IN  SEARCHES 

GEO-OK  -  OVERRIDE  INTERATOMIC  DISTANCE  CHECK 

GNORM=n.n-  FLEPO  EXIT  WHEN  GRADIENT  NORM  BELOW  n.n 

GRADIENTS-  ALL  GRADIENTS  TO  BE  PRINTED 

HCORE  -  PRINT  DETAILS  OF  WORKING  IN  HCORE 

ISOTOPE  -  FORCE  MATRIX  WRITTEN  TO  DISK  (CHANNEL  9  ) 

ITER  -  PRINT  DETAILS  OF  WORKING  IN  ITER 

ITRY-N  -  SET  LIMIT  OF  NUMBER  OF  SCF  ITERATIONS  TO  N. 

LARGE  -  EXPANDED  OUTPUT  TO  BE  PRINTED 
LET  -  DO  NOT  REDUCE  GRADIENTS  IN  FORCE 
LINMIN  -  PRINT  DETAILS  OF  WORKING  IN  LINMIN 
LOCALIZE  -  LOCALIZED  ORBITALS  TO  BE  PRINTED 
LOCMIN  -  PRINT  DETAILS  OF  WORKING  IN  LOCMIN 
MINDO/3  -  THE  MINDO/3  HAMILTONIAN  TO  BE  USED 
MECI  -  PRINT  DETAILS  OF  MECI  CALCULATION 
MOLDAT  -  DETAILS ‘OF  WORKING  IN  MOLDAT  TO  BE  PRINTED 
MULLIK  -  PRINT  THE  MULLIKEN  POPULATION  ANALYSIS 
NLLSQ  -  GRADIENTS  TO  BE  MINIMIZED  USING  NLLSQ. 

NOINTER  -  INTERATOMIC  DISTANCES  NOT  TO  BE  PRINTED 
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NOXYZ  -  CARTESIAN  COORDINATES  NOT  TO  BE  PRINTED 
OLDENS  -  INITIAL  DENSITY  MATRIX  READ  OFF  DISK 
OPEN  -  OPEN-SHELL  RHF  CALCULATION  REQUESTED 
PI  -  RESOLVE  DENSITY  MATRIX  INTO  SIGMA  AND  PI  BONDS 

PL  -  MONITOR  CONVERGANCE  OF  DENSITY  MATRIX  IN  ITER 

PRECISE  -  CRITERIA  TO  BE  INCREASED  BY  100  TIMES 
PULAY  -  PULAY'S  CONVERGER  TO  BE  USED  TO  OBTAIN  A  SCF 
RESTART  -  CALCULATION  RESTARTED 

ROOT-n  -  ROOT  n  TO  BE  OPTIMIZED  IN  A  C.I.  CALCULATION 
ROT-n  -  THE  SYMMETRY  NUMBER  OF  THE  SYSTEM  IS  n. 

SADDLE  -  TRANSITION  STATE  TO  BE  OPTIMIZED 

SCFCRT-.n-  DEFAULT  SCF  CRITERION  REPLACED  BY  THE  VALUE  SUPPLIED 

SHIFT-n  -  A  DAMPING  FACTOR  OF  n  DEFINED 

SINGLET  -  RHF  SINGLET  STATE  REQUIRED 

SPIN  -  FINAL  UHF  SPIN  MATRIX  TO  BE  PRINTED 

SYMMETRY  -  SYMMETRY  CONDITIONS  TO  BE  IMPOSED 

T-n  -  A  TIME  OF  n  SECONDS  REQUESTED 

THERMO  -  A  THERMODYNAMICS  CALCULATION  IS  TO  BE  PERFORMED. 

TIMES  -  TIMES  OF  VARIOUS  STAGES  TO  BE  PRINTED 

TRANS  -  THE  SYSTEM  IS  A  TRANSITION  STATE 

(USED  IN  THERMODYNAMICS  CALCULATION) 

TRIPLET  -  TRIPLET  STATE  REQUIRED 
UHF  -  UNRESTRICTED  HARTREE-FOCK  CALCULATION 
VECTORS  -  FINAL  EIGENVECTORS  TO  BE  PRINTED 
XYZ  -  ALL  GEOMETRIC  OPERATIONS  TO  BE  DONE  IN 
CARTESIAN  COORDINATES. 
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2.3  DEFINITIONS  OF  KEY-WORDS 


1  ELECTRON 

The  final  one-electron  matrix  Is  printed  out.  This  matrix  Is 
composed  of  atomic  orbitals;  the  array  element  between  orbitals  1  and  j 
on  different  atoms  Is  given  by 

H(  1 , j )  -  0.5  x  (betad)  +beta(j))  x  overlapd.j) 

The  matrix  elements  between  orbitals  1  and  j  on  the  same  atom  are 
calculated  from  the  electron-nuclear  attraction  energy,  and  also  from 
the  U<1 )  value,  If  Uj . 

The  one-electron  matrix  Is  unaffected  by  (a)  the  charge,  and  (b) 
the  electron  density.  It  Is  only  a  function  of  the  geometry. 
Abbreviation:  1ELEC. 


The  data  can  be  read  In  and  output,  but  no  actual  calculation  is 
performed  when  this  key-word  Is  used.  This  Is  useful  for  checking 
purposes. 


When  a  single  geometry  Is  to  be  studied,  then  1SCF  should  be  used. 
All  the  key-words  relevant  to  output  can  be  used.  If  the  gradients  are 
to  be  calculated,  then  GRADIENTS  should  be  specified  as  they  are  not 
calculated  by  default. 

If  the  key-word  RESTART  is  also  present,  then  the  geometric 
parameters  which  were  being  optimized  will  be  used  In  the  gradient 
calculation. 

1 SCF  Is  helpful  in  a  learning  situation.  MOPAC  normally  performs 
many  SCF  calculations,  and  In  order  to  minimize  output  when  following 
the  working  of  the  SCF  calculation,  1 SCF  Is  very  useful. 


The  new  AMI  method  Is  to  be  used.  By  default  MNDO  Is  run. 
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BAR- 

In  the  SADDLE  calculation  the  distance  between  the  two  geometries 
Is  steadily  reduced  until  the  transition  state  Is  located.  Sometimes, 
however,  the  user  may  want  to  alter  the  maximum  rate  at  which  the 
distance  between  the  two  geometries  reduces.  BAR  Is  a  ratio,  normally 
0.15,  or  15  percent.  This  represents  a  maximum  rate  of  reduction  of  the 
bar  of  15  percent  per  step.  Alternative  values  that  might  be  considered 
are  BAR-0.05  or  BAR-0.10,  although  other  values  may  be  used.  See  also 
SADDLE. 


BIRADICAL 

NOTE:  BIRADICAL  Is  a  redundant  key-word,  and  represents  a 
particular  configuration  Interaction  calculation.  Experienced  users  of 
MECI  can  duplicate  the  effect  of  the  key-word  BIRADICAL  by  using  the 
MECI  key-words  0PEN<2,2)  and  SINGLET. 

For  molecules  which  are  believed  to  have  blradlcalold  character  the 
option  exists  to  optimize  the  lowest  singlet  energy  state  which  results 
from  the  mixing  of  three  states.  These  states  are,  in  order,  (1)  the 
(mlcro)state  arising  from  a  one  electron  excitation  from  the  HOMO  to  the 
LUMO,  which  Is  combined  with  the  microstate  resulting  from  the 
time-reversal  operator  acting  on  the  parent  microstate,  the  result  being 
a  full  singlet  state;  (2)  the  state  resulting  from  de-excitation  from 
the  formal  LUMO  to  the  HOMO;  and  (3)  the  state  resulting  from  the  single 
electron  In  the  formal  HOMO  being  excited  Into  the  LUMO. 

Microstate  1  Microstate  2  Microstate  3 


Alpha  Beta  Alpha  Beta  Alpha  Beta  Alpha  Beta 


LUMO 


A  configuration  Interaction  calculation  Is  Involved  here.  A  biradical 
calculation  done  without  C.I.  at  the  RHF  level  would  be  meaningless. 
Either  rotational  Invariance  would  be  lost,  as  In  the  D2d  form  of 
ethylene,  or  very  artificial  barriers  to  rotations  would  be  found,  such 
as  In  a  methane  molecule  "orbiting"  a  D2d  ethylene.  In  both  cases  the 
Inclusion  of  limited  configuration  Interaction  corrects  the  error. 
BIRADICAL  should  not  be  used  If  either  the  HOMO  or  LUMO  is  degenerate; 
In  this  case,  the  full  manifold  of  HOMO  x  LUMO  should  be  Included  In  the 
C.I.,  using  MECI  options.  The  user  should  be  aware  of  this  situation. 
When  the  biradical  calculation  Is  performed  correctly,  the  result  Is 
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normally  a  net  stabilization.  However,  If  the  first  singlet  excited 
state  Is  much  higher  In  energy  than  the  closed-shell  ground  state, 
BIRADICAL  can  lead  to  a  destabilization.  Abbreviation:  BIRAD.  See 
also  MECI,  C.I.,  OPEN,  SINGLET. 


BONDS 

The  rotatlonally  Invariant  bond  order  between  all  pairs  of  atoms  Is 
printed.  In  this  context  a  bond  Is  defined  as  the  sum  of  the  squares  of 
the  density  matrix  elements  connecting  any  two  atoms.  For  ethane, 
ethylene,  and  acetylene  the  carbon-carbon  bond  orders  are  roughly  1.00, 
2.00,  and  3.00  respectively.  The  diagonal  terms  are  the  valencies 
calculated  from  the  atomic  terms  only  and  are  defined  as  the  sum  of  the 
bonds  the  atom  makes  with  other  atoms.  In  UHF  and  non-variatlonal ly 
optimized  wavefunctlons  the  calculated  valency  will  be  Incorrect,  the 
degree  of  error  being  proportional  to  the  non-duodempotency  of  the 
density  matrix.  For  an  RHF  wavefunctlon  the  square  of  the  density 
matrix  Is  equal  to  twice  the  density  matrix. 


C. I .=2 

Normally  configuration  Interaction  is  Invoked  If  any  of  the 
key-words  which  Imply  a  C.I.  calculation  are  used,  such  as  BIRADICAL, 
EXCITED  or  TRIPLET.  Note  that  ROOT-  does  not  Imply  a  C.I.  calculation: 
ROOT-  Is  only  used  when  a  C.I.  calculation  Is  done. 

Because  of  storage  limitations,  the  maximum  number  of  M.O.'s 
allowed  In  the  C.I.  Is  2. 

To  form  MOPAC  with  C.I.  use  the  libraries  MOPACLIB  and  MECILIB. 
Associated  key-words:  MECI,  ROOT-,  SINGLET,  DOUBLET,  etc. 


CHARGE- 

When  the  system  being  studied  Is  an  Ion,  the  charge,  n,  on  the  Ion 
can  be  supplied  by  CHARGE-n.  For  cations  n  can  be  1  or  2  or  3,  etc,  for 
anions  -1  or  -2  or  -3,  etc. 

EXAMPLES 


ION 

KEYWORD 

ION 

KEYWORD 

NH4(+) 

CHARGE- 1 

CH3C00(-> 

CHARGE— 1 

C2H5(+) 

CHARGE- 1 

(COO) (-) 

CHARGE— 2 

S04(-) 

CHARGE— 2 

P04<3-> 

CHARGE— 3 

HS04( — ) 

CHARGE— 1 

H2P04(-> 

CHARGE— 1 

KE) -WORDS 


CYCLES- 

In  Bartel's  method  of  gradient  norm  minimization,  NLLSQ,  the 
default  number  of  cycles  (100)  Is  replaced  by  the  number  n  specified  by 
CYCLES-n. 


DCART 

The  cartesian  derivatives  which  are  calculated  in  DCART  for 
variational ly  optimized  systems  are  printed  if  the  key-word  DCART  is 
present.  The  derivatives  are  In  units  of  kcal s/Angstrom,  and  the 
coordinates  are  displacements  in  x,  y,  and  z. 


DEBUG 

Certain  key-words  have  specific  output  control  meanings,  such  as 
FOCK,  VECTORS  and  DENSITY.  If  they  are  used,  only  the  final  arrays  of 
the  relevant  type  are  printed.  If  DEBUG  Is  supplied,  then  all  arrays 
are  printed.  This  is  useful  In  debugging  ITER.  DEBUG  can  also  Increase 
the  amount  of  output  produced  when  a  key-word  Is  used,  e.g.  COMPFG. 


DENOUT 

The  density  matrix  at  the  end  of  the  calculation  is  to  be  output  in 
a  form  suitable  for  input  In  another  job.  If  an  automatic  dump  due  to 
the  time  being  exceeded  occurs  during  the  current  run  then  DENOUT  Is 
Invoked  automatically,  (see  RESTART) 


DENSITY 

At  the  end  of  a  job,  when  the  results  are  being  printed,  the 
density  matrix  Is  also  printed.  For  RHF  the  normal  density  matrix  Is 
printed,  for  UHF  the  addition  of  the  alpha  and  beta  density  matrices  Is 
printed. 

If  density  Is  not  requested,  then  the  diagonal  of  the  density 
matrix,  that  Is,  the  electron  density  on  the  ctomlc  orbitals,  will  be 
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DOUBLET 

When  a  configuration  interaction  calculation  is  done,  all  spin 
states  are  calculated  simultaneously,  either  for  component  of  spin  *  0 
or  1/2.  When  only  doublet  states  are  of  Interest,  then  DOUBLET  can  be 
specified,  and  all  other  spin  states,  while  calculated,  are  Ignored  in 
the  choice  of  root  to  be  used. 

Note  that  while  almost  every  odd-electron  system  will  have  a 
doublet  ground  state,  DOUBLET  should  still  be  specified  if  the  desired 
state  must  be  a  doublet. 

DOUBLET  has  no  meaning  in  a  UHF  calculation. 


DUMP 

Restart  files  are  written  automatically  at  one  hour  cpu  time  intervals 
to  allow  a  long  job  to  be  restarted  if  the  job  is  terminated 
catastrophical ly.  To  change  the  frequency  of  dump,  set  DUMP-nn  to 
request  a  dump  every  nn  seconds.  Alternative  form,  DUMP=nnM  for  a  dump 
every  nn  minutes.  DUMP  only  works  with  geometry  optimization,  gradient 
minimization,  and  FORCE  calculations.  It  does  not  (yet)  work  with  a 
path  or  SADDLE  calculation. 


ECHO 

Data  are  echoed  back  if  ECHO  is  specified.  Only  useful  if  data  are 
suspected  to  be  corrupt. 


ENPART  * 

This  is  a  very  useful  tool  for  analyzing  the  energy  terms  within  a 
system.  The  total  energy,  in  eV,  obtained  by  the  addition  of  the 
electronic  and  nuclear  terms,  is  partitioned  into  mono-  and  bi-centric 
contributions,  and  these  contributions  in  turn  are  divided  into  nuclear 
and  one-  and  two-electron  terms.  ENPART  is  not  included  in  the  minimum 
MOPAC  version.  To  use  ENPART  you  must  link  MOPACLIB  with  ENPARTLIB. 
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ESR 

The  unpaired  spin  density  arising  from  an  odd-electron  system  can 
be  calculated  both  RHF  and  UHF.  In  a  UHF  calculation  the  alpha  and  beta 
M.O.s  have  different  spatial  forms,  so  unpaired  spin  density  can 
naturally  be  present  on  In-plane  hydrogen  atoms  such  as  In  the  phenoxy 
radical . 

In  the  RHF  formalism  a  HECI  calculation  Is  performed.  If  the 
key-words  OPEN  and  C.I.-  are  both  absent  then  only  a  single  state  Is 
calculated.  The  unpaired  spin  density  Is  then  calculated  from  the  state 
function.  In  order  to  have  unpaired  spin  density  on  the  hydrogens  In, 
for  example,  the  phenoxy  radical,  several  states  should  be  mixed. 


EXCITED 

The  state  to  be  calculated  Is  the  first  excited  open-shell  singlet 
state.  If  the  ground  state  Is  a  singlet,  then  the  state  calculated  will 
be  S(1>;  If  the  ground  state  Is  a  triplet,  then  S(2).  This  state  would 
normally  be  the  state  resulting  from  a  one-electron  excitation  from  the 
HOMO  to  the  LUMO.  Exceptions  would  be  If  the  lowest  singlet  state  were 
a  biradical,  In  which  case  the  EXCITED  state  could  be  a  closed  shell. 

The  EXCITED  state  will  be  calculated  from  a  BIRADICAL  calculation 
In  which  the  second  root  of  the  C.I.  matrix  Is  selected.  Note  that  the 
eigenvector  of  the  C.I.  matrix  Is  not  used  In  the  current  formalism. 
Abbreviation:  EXCI. 

NOTE:  EXCITED  Is  a  redundant  key-word,  and  represents  a  particular 
configuration  Interaction  calculation.  Experienced  users  of  MECI  can 
duplicate  the  effect  of  the  '  ey-word  EXCITED  by  using  the  MECI  key-words 
0PEN(2 ,2) ,  SINGLET,  and  ROOT-2. 


FILL- 

The  n'th  M.O.  In  an  RHF  calculation  Is  constrained  to  be  filled. 
It  has  no  effect  on  a  UHF  calculation.  After  the  first  Iteration  (NOTE: 
not  after  the  first  SCF  calculation,  but  after  the  first  Iteration 
within  the  first  SCF  calculation)  the  n’th  M.O.  Is  stored,  and.  If 
occupied,  no  further  action  Is  taken  at  that  time.  If  unoccupied,  then 
the  HOMO  and  the  n'th  M.O.'s  are  swapped  around,  so  that  the  n'th  M.O. 
Is  now  filled.  On  all  subsequent  Iterations  the  M.O.  nearest  In 
character  to  the  stored  M.O.  Is  forced  to  be  filled,  and  the  stored 
M.O.  replaced  by  that  M.O.  This  Is  necessitated  by  the  fact  that  In  a 
reaction  a  particular  M.O.  may  change  Its  character  very  considerably. 
A  useful  procedure  Is  to  run  1SCF  and  DENOUT  first,  In  order  to  Identify 
the  M.O.'s;  the  complete  job  Is  then  run  with  OLDENS  and  FILL-nn,  so 
that  the  eigenvectors  at  the  first  Iteration  are  fully  known.  As  FILL 
Is  known  to  give  difficulty  at  times,  consider  also  using  C.I.-n  and 
ROOT-m. 


-  15  - 


ICtT-WCIKUb 


r aye  c- i v 


FORCE 

A  force-calculation  Is  to  be  run.  The  Hessian,  that  Is  the  matrix 
(In  mlllldynes  per  Angstrom)  of  second  derivatives  of  the  energy  with 
respect  to  displacements  cf  all  pairs  of  atoms  In  x,  y,  and  z,  Is 
calculated.  On  diagonal izatlon  this  gives  the  force  constants  for  the 
molecule.  The  force  matrix,  weighted  for  Isotopic  masses,  Is  then  used 
for  calculating  the  vibrational  frequencies.  The  system  can  be 
characterized  as  a  ground  state  or  a  transition  state  by  the  presence  of 
five  (for  a  linear  system)  or  six  eigenvalues  which  are  very  small  (less 
than  about  30  reciprocal  centimeters).  A  transition  state  is  further 
characterized  by  one,  and  exactly  one,  negative  force  constant. 

A  FORCE  calculation  is  a  prerequisite  for  a  THERMO  calculation. 

Before  a  FORCE  calculation  is  started,  a  check  Is  made  to  ensure 
that  a  stationary  point  is  being  used.  This  check  involves  calculating 
the  gradient  norm  (GNORM)  and  if  it  is  significant,  the  GNORM  will  be 
reduced  using  NLLSQ  (Bartel's  method).  All  internal  coordinates  are 
optimized,  and  any  symmetry  constraints  are  ignored  at  this  point.  An 
implication  of  this  is  that  if  the  specification  of  the  geometry  relies 
on  any  angles  being  exactly  180  or  zero  degrees,  the  calculation  may 
fail. 


The  geometric  definition  supplied  to  FORCE  should  not  rely  on 
angles  or  dihedrals  assuming  exact  values.  (The  test  of  exact  linearity 
Is  sufficiently  slack  that  most  molecules  that  are  linear,  such  as 
acetylene  and  but-2-yne,  should  not  be  stopped  )  See  also  THERMO,  LET, 
TRANS,  ISOTOPE.  Because  of  its  size,  FORCE  calculations  are  not 
Included  in  tne  minimum  MOPAC  PC  version.  To  use  FORCE  you  must  link 
MOPACLIB  witn  FQRCELIB. 


FULSCF 

In  line-searches  the  option  exists  to  require  all  energy 
evaluations  to  be  done  using  full  SCF  calculations.  Normally  full  SCF 
calculations  are  not  carried  out  during  a  line  search  as  the  density 
matrix  is  normally  not  changing  very  fast.  The  only  Important  exception 
is  In  non-variational ly  optimized  wavefuncttons ,  such  as  occur  In 
half-electron  or  C.I.  calculations. 


GE0-0K 

Normally  the  program  will  stop  with  a  warning  message  if  two  atoms 
are  within  0.8  Angstroms  of  each  other,  or  (more  rarely)  the  D.F.P. 
routine  has  difficulty  optimizing  the  geometry.  GEO-OK  will  over-ride 
the  job  termination  sequence,  and  allow  the  calculation  to  proceed.  In 
practice  most  jobs  that  terminate  due  to  these  checks  contain  errors  In 
data,  so  caution  should  be  exercised  if  GEO-OK  is  used.  An  Important 
exception  to  this  warning  is  if  the  system  contains,  or  may  give  rise 
to,  a  Hydrogen  molecule.  GEO-OK  will  override  other  geometric  safety 
checks  such  as  the  unstable  gradient  in  a  geometry  optimization 
prevent. ng  reliable  optimization. 
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See  also  the  message  "GRADIENTS  OF  OLD  GEOMETRY,  GNORM=  nn.nnnn" 


GNORM- 

The  D.F.P.  geometry  optimization  termination  criteria  can  be 
over-ridden  by  specifying  a  gradient  norm  requirement.  For  example, 
GN0RM-20  would  allow  the  D.F.P.  to  exit  as  soon  as  the  gradient  norm 
dropped  below  20.0,  the  default  being  1.0.  A  GNORM-O.Ol  could  be  used 
to  refine  a  geometry  beyond  the  normal  limits.  WARNING:  If  a  very 
small  value  Is  chosen,  the  D.F.P.  procedure  may  not  terminate  In  a 
reasonable  time.  A  reasonable  lower  bound  for  GNORM  Is  0.1. 


GRADIENTS 

In  a  1SCF  calculation  gradients  are  not  calculated  by  default:  In 
non-varlatlonal ly  optimized  systems  this  would  take  an  excessive  time. 
GRADIENTS  allows  the  gradients  to  be  calculated.  All  gradients  are  then 
calculated,  whether  marked  for  calculation  or  not,  and  printed.  An 
exception  Is  when  the  1SCF  was  used  In  conjunction  with  the  key-word 
RESTART,  In  which  case  only  the  coordinates  being  optimized  would  have 
their  gradients  printed.  Abbreviation:  GRAD. 


ISOTOPE 

The  FORCE  matrix  Is  very  time-consuming  to  generate,  and  In 
Isotopic  substitution  studies  several  vibrational  calculations  may  be 
needed.  To  allow  the  frequencies  to  be  calculated  from  the  (constant) 
force  matrix,  ISOTOPE  Is  used.  When  a  FORCE  calculation  is  completed, 
ISOTOPE  will  cause  the  force  matrix  to  be  stored,  regardless  of  whether 
or  not  any  Intervening  restarts  have  been  made.  To  re-calculate  the 
frequencies,  etc.,  starting  at  the  end  of  the  force  matrix  calculation 
specify  RESTART. 

The  two  key-words  RESTART  and  ISOTOPE  can  be  used  together.  For 
example,  If  a  normal  FORCE  calculation  runs  for  a  long  time,  the  user 
may  want  to  divide  It  up  Into  stages  and  save  the  final  force  matrix. 
Once  ISOTOPE  has  been  used.  It  does  not  need  to  be  used  on  subsequent 
RESTART  runs. 


ITRY=NN 

The  default  maximum  number  of  SCF  Iterations  is  200.  When  this 
limit  presents  difficulty,  ITRY=nn  can  be  used  to  re-deflne  it.  For 
example,  If  ITRY-400  Is  used,  the  maximum  number  of  Iterations  will  be 
set  to  400.  ITRi  should  normally  not  be  changed  until  all  other  means 
of  obtaining  a  SCF  have  been  exhausted,  e.g.  SHIFT  and  PULAY  etc. 
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LET 

Before  the  Hessian  matrix  Is  calculated  in  a  FORCE  calculation  the 
geometry  will  be  refined  by  a  gradient  minimization  routine  if  the 
gradient  norm  is  significant.  If  the  user  does  not  want  the  refinement 
to  be  carried  out,  then  LET  is  provided  to  let  the  calculation  proceed. 


LOCALIZE 

The  occupied  eigenvectors  are  transformed  into  a  localized  set  of 
M.O.'s  by  a  series  of  2  by  2  rotations  which  maximize  < ps i * *4 > .  The 
value  of  l/<psi**4>  is  a  direct  measure  of  the  number  of  centers 
Involved  in  the  M.O.,  thus  for  H2  the  value  of  1 / <ps i **4>  is  2.0,  for  a 
three-center  bond  is  3.0,  and  a  lone  pair  would  be  1.0.  Higher 
degeneracies  than  allowed  by  point  group  theory  are  readily  obtained. 
For  example,  benzene  would  give  rise  to  a  6-fold  degenerate  C-H  bond,  a 
6-fold  degenerate  C-C  sigma  bond  and  a  three-fold  degenerate  C-C  pi 
bond.  In  principle,  there  is  no  single  step  method  to  unambiguously 
obtain  the  most  localized  set  of  M.O.'s  in  systems  where  several 
canonical  structures  are  possible,  just  as  no  simple  method  exists  for 
finding  the  most  stable  conformer  of  some  large  compound.  However,  the 
localized  bonds  generated  will  normally  be  quite  acceptable  for  routine 
applications.  Abbreviation:  LOCAL. 


MECI 

At  the  end  of  the  calculation  details  of  the  Multi  Electron 
Interaction  Calculation  are  printed  if  MECI  is  specified.  The  state 
vectors  can  be  printed  by  specifying  VECTORS.  The  MECI  calculation  is 
either  invoked  automatically,  or  explicitly  Invoked  by  the  use  of  the 
C.I.=n  keyword. 
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MICROS-n 

The  microstates  used  by  MECI  are  normally  generated  by  use  of  a 
permutation  operator.  When  Individually  defined  microstates  are 
desired,  then  MICROS-n  can  be  used,  where  n  defines  the  number  of 
microstates  to  be  read  In. 

Format  for  Microstates 

After  the  geometry  data  plus  any  symmetry  data  are  read  In,  data 
defining  each  microstate  Is  read  in,  using  format  2011,  one  microstate 
per  line. 

For  a  system  with  n  M.O.s  In  the  C.I.  (use  OPEN-(nl.n)  or  C.I.-n 
to  do  this),  the  populations  of  the  n  alpha  M.O.s  are  defined,  followed 
by  the  n  beta  M.O.s.  Allowed  occupancies  are  zero  and  one.  For  n-2  the 
closed-shell  ground  state  would  be  defined  as  1010,  meaning  one  electron 
In  the  alpha  H.O.M.O.  and  one  electron  In  the  beta  H.O.M.O. 

Users  are  warned  that  they  are  responsible  for  completing  any  spin 
manifolds.  Thus  while  the  state  1100  Is  a  triplet  state  with  component 
of  spin  -  1,  the  state  1001,  while  having  a  component  of  spin  =  0  Is 
neither  a  singlet  nor  a  triplet.  In  order  to  complete  the  spin  manifold 
the  microstate  0110  must  also  be  Included. 

If  a  manifold  of  spin  states  Is  not  complete,  then  the  eigenstates 
of  the  spin  operator  will  not  be  quantized.  When  4  or  fewer  microstates 
are  supplied,  this  Is  the  only  cause  of  loss  of  spin  quantization. 

There  are  two  other  limitations  on  possible  microstates.  First, 
the  number  of  electrons  In  every  microstate  should  be  the  same.  If  they 
differ,  a  warning  message  will  be  printed,  and  the  calculation  continued 
(but  the  results  will  almost  certainly  be  nonsense).  Second,  the 
component  of  spin  for  every  microstate  must  be  the  same,  except  for 
teaching  purposes.  Two  microstates  of  different  components  of  spin  will 
have  a  zero  matrix  element  connecting  them.  No  warning  will  be  given  as 
this  Is  a  reasonable  operation  In  a  teaching  situation.  For  example,  If 
all  states  arising  from  two  electrons  In  two  levels  are  to  be 
calculated,  say  for  teaching  Russel-Saunders  coupling,  then  the 
following  microstates  would  be  used: 
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Microstate  No.  of  alpha,  beta  electrons  Ms  State 

1100  2  0  1  Triplet 

1010  1  1  0  Singlet 

1001  1  1  0  Mixed 

0110  1  1  0  Mixed 

0101  1  1  0  Singlet 

0011  02  -1  Triplet 

Constraints  on  the  space  manifold  are  just  as  rigorous,  but  much 
easier  to  satisfy.  If  the  energy  levels  are  degenerate,  then  all 
components  of  a  manifold  of  degenerate  M.O.s  should  be  either  Included 
or  excluded.  If  only  some,  but  not  all,  components  are  used,  the 
required  degeneracy  of  the  states  will  be  missing. 

The  format  is  defined  as  2011  so  that  spaces  can  be  used  for  empty 
M.O.s. 


MINDO/3 

The  default  Hamiltonian  within  M0PAC  is  MNDO,  with  the  alternatives 
of  AMI  and  MINDO/3.  To  use  the  MINDO/3  Hamiltonian  the  key-word  MINDO/3 
should  be  used.  Acceptable  alternatives  to  the  key-word  MINDO/3  are 
MINDO  and  MIND03. 


MULLIK 

A  full  Mulliken  Population  analysis  is  to  be  done  on  the  final  RHF 
wavefunction.  This  Involves  the  following  steps: 

(1)  The  eigenvector  matrix  is  divided  by  the  square  root 
of  the  overlap  matrix,  S. 

(2)  The  Coul son-type  density  matrix,  P,  is  formed. 

(3)  The  overlap  population  is  formed  from  P( 1 , j )*S( 1 , j ) . 

(4)  Half  the  off-diagonals  are  added  onto  the  diagonals. 

MULLIK  is  not  supported  when  using  FORCE  or  NLLSQ-  To  use  MULLIK  you 
must  link  MOPACLIB  with  ENPARTLIB. 
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NLLSQ 

The  gradient  norm  Is  to  be  minimized  by  Bartel's  method.  This  Is  a 
Non-Linear  Least  Squares  gradient  minimization  routine.  Gradient 
minimization  will  locate  one  of  three  possible  points: 

(a)  A  minimum  In  the  energy  surface.  The  gradient  norm  will  go  to 
zero,  and  the  lowest  five  or  six  eigenvalues  resulting  from  a  FORCE 
calculation  will  be  approximately  zero. 

<b>  A  transition  state.  The  gradient  norm  will  vanish,  as  In  (a), 
but  In  this  case  the  system  Is  characterized  by  one,  and  only  one, 
negative  force  constant. 

<c)  A  local  minimum  In  the  gradient  norm  space.  In  this  (normally 
unwanted)  case  the  gradient  norm  is  minimized,  but  does  not  go  to  zero. 
A  FORCE  calculation  will  not  give  the  five  or  six  zero  eigenvalues 
characteristic  of  a  stationary  point.  While  normally  undesirable,  this 
Is  sometimes  the  only  way  to  obtain  a  geometry.  For  Instance,  If  a 
system  Is  formed  which  cannot  be  characterized  as  an  Intermediate,  and 
at  the  same  time  Is  not  a  transition  state,  but  nonetheless  has  some 
chemical  significance,  then  that  state  can  be  refined  using  NLLSQ. 
Because  of  Its  size,  NLLSQ  Is  not  supported  when  using  the  minimal  MOPAC 
PC  Version.  To  use  NLLSQ  you  must  link  with  NLLSQLIB.  When  NLLSQ  Is 
Included,  the  keywords  BONDS,  MULLIK,  and  LOCAL  are  not  supported. 


OLDENS 

A  density  matrix  produced  by  an  earlier  run  of  MOPAC  Is  to  be  used 
to  start  the  current  calculation.  This  can  be  used  In  attempts  to 
obtain  a  SCF  when  a  previous  calculation  ended  successfully  but  a 
subsequent  run  failed  to  go  SCF. 


OPENCnl ,n2) 

The  M.O.  occupancy  during  the  SCF  calculation  can  be  defined  In 
terms  of  doubly  occupied,  empty,  and  fractionally  occupied  M.O.s.  The 
fractionally  occupied  M.O.s  are  defined  by  0PEN(nl,n2),  where  nl  * 
number  of  electrons  In  the  open-shell  manifold,  and  n2  =  number  of 
open-shell  M.O.s;  nl  must  be  In  the  range  0  to  2.  OPEN(l.l)  will  be 
assumed  for  odd-electron  systems  unless  an  OPEN  keyword  is  used.  Errors 
Introduced  by  use  of  fractional  occupancy  are  automatically  corrected  In 
a  MECI  calculation  when  0PEN(nl,n2)  is  used. 
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PI 

The  normal  density  matrix  Is  composed  of  atomic  orbitals,  that  Is 
s,  px,  py  and  pz.  PI  allows  the  user  to  see  how  each  atom-atom 
Interaction  Is  split  Into  sigma  and  pi  bonds.  The  resulting  "density 
matrix"  Is  composed  of  the  following  basis-functions:-  s-sigma,  p-slgma, 
p-pl ,  d-slgma,  d-pl,  d-dell.  The  on-diagonal  terms  give  the 
hybridization  state,  so  that  an  sp2  hybridized  system  would  be 
represented  as  s-sigma  1.0,  p-slgma  2.0,  p-pl  1.0 


k 


: 


PRECISE 


The  criteria  for  terminating  all  optimizations,  electronic  and 
geometric,  are  to  be  increased  by  a  factor,  normally,  100.  This  can  be 
used  where  more  precise  results  are  wanted.  If  the  results  are  going  to 
be  used  In  a  FORCE  calculation,  where  the  geometry  needs  to  be  known 
quite  precisely,  then  PRECISE  is  recommended;  for  small  systems  the 
extra  cost  In  CPU  time  Is  minimal. 


PULAY 


The  default  converger 
Pulay’s  procedure  as  soon 
A  considerable  improvement 
If  a  large  number  of 
calculation  using  1SCF  and 
Its  own,  and  if  a  saving 
the  full  calculation. 


In  the  SCF  calculation  is  to  be  replaced  by 
as  the  density  matrix  Is  sufficiently  stable, 
in  speed  can  be  achieved  by  the  use  of  PULAY. 
SCF  calculations  are  envisaged,  a  sample 
PULAY  should  be  compared  with  using  1SCF  on 
In  time  results,  then  PULAY  should  be  used  in 


RESTART 


When  a  job  has  been  stopped,  for  whatever  reason,  and  Intermediate 
results  have  been  stored,  then  the  calculation  can  be  restarted  at  the 
point  where  It  stopped  by  specifying  RESTART.  The  most  common  cause  of 
a  job  stopping  before  completion  Is  its  exceeding  the  time  allocated.  A 
saddle-point  calculation  has  no  restart,  but  the  output  file  contains 
Information  which  can  easily  be  used  to  start  the  calculation  from  a 
point  near  to  where  It  stopped. 


It  Is  not  necessary  to  change  the  geometric  data  to  reflect  the  new 
geometry,  as  a  result  the  geometry  printed  at  the  start  of  a  restarted 
job  will  be  that  of  the  original  data,  not  that  of  the  restarted  file. 


A  convenient  way  to  monitor  a  long  run  Is  to  specify  1 SCF  and 
RESTART;  this  will  give  a  normal  output  file  at  very  little  cost.  NOTE 
1:  In  the  FORCE  calculation  two  restarts  are  possible.  These  are  (a)  a 
restart  In  FLEP0  If  the  geometry  was  not  optimized  fully  before  FORCE 
was  called,  and  (b)  the  normal  restart  in  the  construction  of  the  force 
matrix.  If  the  restart-  Is  in  FLEP0  within  FORCE  then  the  key-word  FORCE 
should  be  deleted,  and  the  key-word  RESTART  used  on  its  own. 
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ROOT-n 

The  n'th  root  of  a  C.I.  calculation  Is  to  be  used  In  the 
calculation.  If  a  key-word  specifying  the  spin-state  Is  also  present, 
e.g.  SINGLET  or  TRIPLET,  then  the  n’th  root  of  that  state  will  be 
selected.  Thus  ROOT-3  and  SINGLET  will  select  the  third  singlet  root. 
If  ROOT-3  Is  used  on  Its  own,  then  the  third  root  will  be  used,  which 
may  be  a  triplet,  the  third  singlet,  or  the  second  singlet  (the  second 
root  might  be  a  triplet).  In  normal  use,  this  key-word  would  not  be 
used.  It  Is  retained  for  educational  and  research  purposes.  Unusual 
care  should  be  exercised  when  ROOT-  Is  specified. 


ROT-n 

In  the  calculation  of  the  rotational  contributions  to  the 
thermodynamic  quantities  the  symmetry  number  of  the  molecule  must  be 
supplied.  The  symmetry  number  of  a  point  group  Is  the  number  of 
equivalent  positions  attainable  by  pure  rotations.  No  reflections  or 
Improper  rotations  are  allowed.  This  number  cannot  be  assumed  by 
default,  and  may  be  affected  by  subtle  modifications  to  the  molecule, 
such  as  isotopic  substitution.  A  list  of  the  most  Important  symmetry 
numbers  follows: 


TABLE  OF  SYMMETRY  NUMBERS 


Cl 

Cl  CS 

1 

D2 

D2D 

C2 

C2V 

C2H 

2 

D3 

030 

C3 

C3V 

C3H 

3 

D4 

D4D 

C4 

C4V 

C4H 

4 

D6 

D6D 

C6 

C6V 

C6H 

6 

S6 

D2H 

4 

C(INF>V 

1 

03H 

6 

D(INF>H 

2 

D4H 

8 

T  TD 

12 

D6H 

12 

OH 

24 

3 


SADDLE 


The  transition  state  In  a  simple  chemical  reaction  Is  to  be 
optimized.  Extra  data  are  required.  After  the  first  geometry, 
specifying  the  reactants,  and  any  symmetry  functions  have  been  defined, 
the  second  geometry,  specifying  the  products,  Is  defined,  using  the  same 
format  as  that  of  the  first  geometry. 


SADDLE  often  falls  to  work  successfully.  Frequently  this  Is  due  to 
dihedral  angles  being  different  by  about  360  degrees  rather  than  zero 
degrees.  As  the  choice  of  dihedral  can  be  difficult,  users  should 
consider  running  this  calculation  with  the  key-word  XYZ.  There  is 
normally  no  ambiguity  In  the  definition  of  cartesian  coordinates.  See 
also  BAR-. 


Many  of  the  bugs  In  SADDLE  have  been  removed  In  this  version.  Use 
of  the  XYZ  option  Is  strongly  recommended.  SADDLE  Is  not  supported  In 
the  minimum  MOPAC  PC  Version.  To  use  SADDLE  link  MOPACLIB  with 
SADDLE. LIB. 
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SCFCRT- 

The  default  SCF  criterion  is  to  be  replaced  by  that  defined  by 
SCFCRT-.  The  SCF  criterion  can  be  varied  from  about  0.0001  to 
0.0000000001.  To  find  a  suitable  value  1SCF  and  various  values  of 
SCFCRT-n.nnn  should  be  used;  a  SCFCRT  which  allows  evaluation  of  the 
heat  of  formation  to  an  acceptable  precision  can  thus  be  found  rapidly. 
An  overly  tight  criterion  can  lead  to  failure  to  achieve  a  SCF,  and 
consequent  failure  of  the  run. 


SHIFT- 

In  an  attempt  to  obtain  an  SCF  by  damping  oscillations  which  are 
slowing  down  the  convergence  or  preventing  a  SCF  being  achieved,  SHIFT 
can  be  used.  The  principle  Is  that  If  the  virtual  M.O.'s  are  raised  In 
energy  relative  to  the  occupied  set,  then  the  polarizability  of  the 
occupied  M.O. 1 s  wl 1 1  decrease,  the  oscillations  being  attributed  to 
autoregeneratlve  charge  fluctuations.  A  SHIFT  of  20  will  raise  the 
virtual  M.O.'s  by  20  eV  above  their  correct  value.  The  disadvantage  of 
SHIFT  Is  that  a  large  value  can  lead  to  excessive  damping,  and  thus 
prevent  an  SCF  being  generated.  As  some  virtual  M.O.'s  are  used  in 
non-variational ly  optimized  calculations  SHIFT  Is  automatically  annulled 
at  the  end  of  the  SCF  in  these  circumstances.  All  effects  of  SHIFT  are 
removed  before  the  results  are  printed.  SHIFT  should  be  used  with  care 
In  that  the  use  of  any  value  for  SHIFT  will  prevent  the  combined  package 
of  convergers  (SHIFT  and  PULAY)  from  automatically  being  used  in  event 
that  the  system  falls  to  go  SCF  In  < I  TRY- 10)  Iterations. 


SINGLET 

When  a  configuration  Interaction  calculation  Is  done,  all  spin 
states  are  calculated  simultaneously,  either  for  component  of  spin  -  0 
or  1/2.  When  only  singlet  states  are  of  interest,  then  SINGLET  can  be 
specified,  and  all  other  spin  states,  while  calculated,  are  Ignored  in 
the  choice  of  root  to  be  used. 

Note  that  while  almost  every  even-electron  system  will  have  a 
singlet  ground  state,  SINGLET  should  still  be  specified  if  the  desired 
state  must  be  a  singlet. 

SINGLET  has  no  meaning  In  a  UHF  calculation,  but  see  also  TRIPLET. 


SPIN 

The  spin  matrix,  defined  as  the  difference  between  the  alpha  and 
beta  density  matrices,  is  to  be  printed.  If  the  system  has  a 
closed-shell  ground  state,  e.g.  methane  run  UHF,  the  spin  matrix  will 
be  null . 

If  SPIN  is  not  requested  in  a  UHF  calculation,  then  the  diagonal  of 
the  spin  matrix,  that  Is  the  spin  density  on  the  atomic  orbitals,  will 
be  printed. 
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SYMMETRY 

Symmetry  data  defining  related  bond  lengths,  angles  and  dihedrals 
can  be  Included  by  supplying  additional  data  after  the  geometry  has  been 
entered.  If  there  are  any  other  data,  such  as  values  for  the  reaction 
coordinates,  or  a  second  geometry,  as  required  by  SADDLE,  then  it  would 
follow  the  symmetry  data.  Symmetry  data  are  terminated  by  one  blank 
line.  For  non-varlatlonal ly  optimized  systems  symmetry  constraints  can 
save  a  lot  of  time  because  many  derivatives  do  not  need  to  be 
calculated.  At  the  same  time,  there  Is  a  risk  that  the  geometry  may  be 
wrongly  specified,  e.g.  If  methane  radical  cation  Is  defined  as  being 
tetrahedral,  no  Indication  that  this  Is  faulty  will  be  given  until  a 
FORCE  calculation  Is  run.  (This  system  undergoes  spontaneous 

Jahn-Teller  distortion.) 

Usually  a  lower  heat  of  formation  can  be  obtained  when  SYMMETRY  Is 
specified.  To  see  why,  consider  the  geometry  of  benzene.  If  no 
assumptions  are  made  regarding  the  geometry,  then  all  the  C-C  bond 
lengths  will  be  very  slightly  different,  and  the  angles  will  be  almost, 
but  not  quite  120  degrees.  Fixing  all  angles  at  120  degrees,  dihedrals 
at  180  or  0  degrees,  and  only  optimizing  one  C-C  and  one  C-H  bond-length 
will  result  In  a  2-D  optimization,  and  exact  D6h  symmetry.  Any 
deformation  from  this  symmetry  must  Involve  error,  so  by  Imposing 
symmetry  some  error  Is  removed. 

The  layout  of  the  symmetry  data  Is: 
xdeflnlng  atom>  <symmetry  relation)  cdeflned  atom>  xdeflned  atom>,... 
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For  example,  ethane,  with  three  Independent  variables,  can  be 
defined  as 


SYMMETRY 
ETHANE,  D3D 


c 

0.000000  0 

0.000000  0 

0.000000  0 

0 

0 

0 

c 

1.528853  1 

0.000000  0 

0.000000  0 

1 

0 

0 

H 

1.105161  1 

110.240079  1 

0.000000  0 

2 

1 

0 

H 

1.105161  0 

110.240079  0 

120.000000  0 

2 

1 

3 

H 

1.105161  0 

110.240079  0 

240.000000  0 

2 

1 

3 

H 

1.105161  0 

110.240079  0 

60.000000  0 

1 

2 

3 

H 

1.105161  0 

110.240079  0 

180.000000  0 

1 

2 

3 

H 

1.105161  0 

110.240079  0 

300.000000  0 

1 

2 

3 

0 

0.000000  0 

0.000000  0 

0.000000  0 

0 

0 

0 

3, 

1,  4, 

5,  6, 

7,  8. 

3, 

2.  4, 

5,  6, 

7,  8, 

Here  atom  3,  a  hydrogen,  is  used  to  define  the  bond  length 
(function  1)  of  atoms  4, 5, 6, 7  and  8;  similarly,  its  angle  (function  2) 
is  used  to  define  the  bond-angle  of  atoms  4, 5, 6, 7  and  8.  The  other 
angles  are  point-group  symmetry  defined  as  a  multiple  of  60  degrees. 

Spaces,  tabs  or  commas  can  be  used  to  separate  data.  Note  that 
only  three  parameters  are  marked  to  be  optimized.  To  end  symmetry  data 
use  a  blank  line. 
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The  full  list  of  available  symmetry  relations  Is  as  follows: 


SYMMETRY  FUNCTIONS 


1  BOND  LENGTH  IS  SET  EQUAL  TO  THE  REFERENCE  BOND  LENGTH 

2  BOND  ANGLE  IS  SET  EQUAL  TO  THE  REFERENCE  BOND  ANGLE 

3  DIHEDRAL  ANGLE  IS  SET  EQUAL  TO  THE  REFERENCE  DIHEDRAL  ANGLE 

4  DIHEDRAL  ANGLE  VARIES  AS  90  DEGREES  -  REFERENCE  DIHEDRAL 

5  DIHEDRAL  ANGLE  VARIES  AS  90  DEGREES  +  REFERENCE  DIHEDRAL 

6  DIHEDRAL  ANGLE  VARIES  AS  120  DEGREES  -  REFERENCE  DIHEDRAL 

7  DIHEDRAL  ANGLE  VARIES  AS  120  DEGREES  +  REFERENCE  DIHEDRAL 

8  DIHEDRAL  ANGLE  VARIES  AS  180  DEGREES  -  REFERENCE  DIHEDRAL 

9  DIHEDRAL  ANGLE  VARIES  AS  180  DEGREES  +  REFERENCE  DIHEDRAL 

10  DIHEDRAL  ANGLE  VARIES  AS  240  DEGREES  -  REFERENCE  DIHEDRAL 

11  DIHEDRAL  ANGLE  VARIES  AS  240  DEGREES  +  REFERENCE  DIHEDRAL 

12  DIHEDRAL  ANGLE  VARIES  AS  270  DEGREES  -  REFERENCE  DIHEDRAL 

13  DIHEDRAL  ANGLE  VARIES  AS  270  DEGREES  +  REFERENCE  DIHEDRAL 

14  DIHEDRAL  ANGLE  VARIES  AS  THE  NEGATIVE  OF  THE  REFERENCE 
DIHEDRAL 

15  BOND  LENGTH  VARIES  AS  HALF  THE  REFERENCE  BOND  LENGTH 

16  BOND  ANGLE  VARIES  AS  HALF  THE  REFERENCE  BOND  ANGLE 

17  BOND  ANGLE  VARIES  AS  180  DEGREES  -  REFERENCE  BOND  ANGLE 

18  BOND  LENGTH  IS  A  MULTIPLE  OF  REFERENCE  BOND-LENGTH 

Function  18  Is  Intended  for  use  In  polymers.  In  which  the 
translation  vector  may  be  a  multiple  of  some  bond-length.  1,2,3  and  14 
are  most  commonly  used.  Abbreviation:  SYM. 

SYMMETRY  Is  not  available  for  use  with  cartesian  coordinates. 


T- 

Thls  Is  a  facility  to  allow  the  program  to  shut  down  In  an  orderly 
manner  on  computers  with  execution  time  C.P.U.  limits. 

The  total  C.P.U.  time  allowed  for  the  current  job  Is  limited  to 
nn.nn  seconds;  by  default  this  Is  one  hour,  l.e.,  3600  seconds.  If  the 
next  cycle  of  the  calculation  cannot  be  completed  without  running  a  risk 
of  exceeding  the  assigned  time  the  calculation  will  write  a  restart  file 
and  then  stop.  The  safety  margin  Is  100  percent;  that  Is,  to  do  another 
cycle,  enough  time  to  do  at  least  two  full  cycles  must  remain. 

An  alternative  specification  Is  T*nn.nnM,  this  allows  time  In 
minutes  to  be  defined,  thus  T-60M  will  define  3600  seconds,  or  60 
minutes. 
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THERMO 

The  thermodynamic  quantities  Internal  energy,  heat  capacity, 
partition  function,  and  entropy  can  be  calculated  for  translation, 
rotation  and  vibrational  degrees  of  freedom  for  a  single  temperature,  or 
a  range  of  temperatures.  Special  situations  such  as  linear  systems  and 
transition  states  are  accommodated.  The  approximations  used  In  the 
THERMO  calculation  are  invalid  below  100K,  and  checking  of  the  lower 
bound  of  the  temperature  range  Is  done  to  prevent  temperatures  of  less 
than  100K  being  used. 

Another  limitation,  for  which  no  checking  Is  done.  Is  that  there 
should  be  no  Internal  rotations.  If  any  exist,  they  will  not  be 
recognized  as  such,  and  the  calculated  quantities  will  be  too  low  as  a 
result. 

If  THERMO  Is  specified  on  Its  own,  then  the  default  values  of  the 
temperature  range  are  assumed.  This  starts  at  200K  and  Increases  in 
steps  of  10  degrees  to  400K.  Three  options  exist  for  overriding  the 
default  temperature  range.  These  are: 


THERMO(nnn) 

The  thermodynamic  quantities  for  a  200  degree  range  of 
temperatures,  starting  at  nnnK  and  with  an  interval  of  10  degrees  are  to 
be  calculated. 


THERMO(nnn.mmm) 

The  thermodynamic  quantities  for  the  temperature  range  limited  by  a 
lower  bound  of  nnn  Kelvin  and  an  upper  bound  of  mmm  Kelvin,  the  step 
size  being  calculated  In  order  to  give  approximately  20  points,  and  a 
reasonable  value  for  the  step.  The  size  of  the  step  In  Kelvin  degrees 
will  be  1,  2,  or  5,  or  a  power  of  10  times  these  numbers. 


THERMO(nnn,mmm, 1 1 1 > 

As  for  THERMO(nnn,mmm)  only  now  the  user  can  explicitly  define  the 
step  size.  The  step  size  cannot  be  less  than  IK. 
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TRANS 

The  Imaginary  frequency  due  to  the  reaction  vector  In  a  transition 
state  calculation  must  not  be  Included  In  the  thermochemical 
calculation.  The  number  of  genuine  vibrations  considered  can  be: 

3N-5  for  a  linear  ground  state  system, 

3N-6  for  a  non-linear  ground  state  system,  or 

3N-6  for  a  linear  transition-state  complex, 

3N-7  for  a  non-linear  transition-state  complex. 

This  key-word  must  be  used  In  conjunction  with  THERMO  If  a 
transition  state  Is  being  calculated. 


TRANS-n 

The  facility  exists  to  allow  the  THERMO  calculation  to  handle 
systems  with  Internal  rotations.  TRANS-n  will  remove  the  n  lowest 
vibrations.  Note  that  TRANS-1  Is  equivalent  to  TRANS  on  Its  own.  For 
xylene,  for  example,  TRANS-2  would  be  suitable. 


TRIPLET 

The  triplet  state  Is  defined.  If  the  system  has  an  odd  number  of 
electrons,  an  error  message  will  be  printed. 

UHF  Interpretation. 

The  number  of  alpha  electrons  exceeds  that  of  the  beta  electrons  by 
2.  If  TRIPLET  Is  not  specified,  then  the  numbers  of  alpha  and  beta 
electrons  are  set  equal.  This  does  not  necessarily  correspond  to  a 
singlet. 

RHF  Interpretation. 

A  RHF  MECI  calculation  Is  performed  to  calculate  the  triplet  state. 
If  no  other  C.I.  key-words  are  used,  then  only  one  state  Is  calculated 
by  default.  The  occupancy  of  the  M.O.s  In  the  SCF  calculation  Is 
defined  as  < . . .2,1 , 1 ,0, . . ) ,  that  Is,  one  electron  Is  put  In  each  of  the 
two  highest  occupied  M.O.s. 

See  key-words  C.I.«n  and  0PEN(nl,n2). 


UHF 


The  unrestricted  Hartree-Fock  Hamiltonian  Is  to  be  used. 
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VECTORS 

The  eigenvectors  are  to  be  printed.  In  UHF  calculations  both  alpha 
and  beta  eigenvectors  are  printed;  In  all  cases  the  full  set,  occupied 
and  virtual,  are  output.  The  eigenvectors  are  normalized  to  unity,  that 
Is  the  sum  of  the  squares  of  the  coefficients  Is  exactly  one.  If  DEBUG 
Is  specified,  then  ALL  eigenvectors  on  every  Iteration  of  every  SCF 
calculation  will  be  printed.  This  Is  useful  In  a  learning  context,  but 
would  normally  be  very  undesirable. 


The  SADDLE  calculation  quite  often  falls  due  to  faulty  definition 
of  the  second  geometry  because  the  dihedrals  give  a  lot  of  difficulty. 
To  make  this  option  easier  to  use,  XYZ  was  developed.  Using  XY2  the 
calculation  runs  entirely  In  cartesian  coordinates,  this  eliminating  the 
problems  associated  with  dihedrals.  The  connectivity  of  the  two  systems 
can  be  different,  but  the  numbering  must  be  the  same.  Dummy  atoms  can 
be  used;  these  will  be  removed  at  the  start  of  the  run.  A  new  numbering 
system  will  be  generated  by  the  program,  when  necessary. 


CHAPTER  3 


GEOMETRY  SPECIFICATION 


FORMAT:  The  geometry  Is  read  in  using  essentially  "Free-Format"  of 
FORTRAN-77.  In  fact,  a  character  Input  Is  used  In  order  to  accommodate 

the  chemical  symbols,  but  the  numeric  data  can  be  regarded  as 

"free-format" .  This  means  that  Integers  and  real  numbers  can  be 
Interspersed,  numbers  can  be  separated  by  one  or  more  spaces,  a  tab 
and/or  by  one  comma.  If  a  number  Is  not  specified.  Its  value  Is  set  to 
zero. 

The  geometry  can  be  defined  In  terms  of  either  Internal  or 
cartesian  coordinates. 

INTERNAL  COORDINATE  DEFINITION 

For  any  one  atom  (1)  this  consists  of  an  Interatomic  distance  In 
Angstroms  from  an  already-defined  atom  (j),  an  Interatomic  angle  In 

degrees  between  atom  1  and  j  and  an  already  defined  k,  (k  and  j  must  be 

different  atoms),  and  finally  a  torsional  angle  in  degrees  between  atom 
1,  j,  k,  and  an  already  defined  atom  1  (1  cannot  be  the  same  as  k  or  j) 

Exceptions: 

1.  Atom  1  has  no  coordinates  at  all:  this  Is  the  origin. 

2.  Atom  2  must  be  connected  to  atom  1  by  an  Interatomic  distance 
only. 

3.  Atom  3  can  be  connected  to  atom  1  or  2,  and  must  make  an  angle 

with  atom  2  or  1  (thus  -  3-2-1  or  3-1-2);  no  dihedral  is 

possible  for  atom  3.  By  default,  atom  3  Is  connected  to  atom 

2. 
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3.1  CONSTRAINTS 


1.  Interatomic  distances  must  be  greater  than  zero.  Zero 
Angstroms  Is  not  accjptable.  The  only  exception  Is  if  the 
parameter  Is  symmetry-related  to  another  atom,  and  Is  the 
dependent  function. 

2.  Angles  must  be  positive.  This  constraint  Is  for  the  benefit  of 
the  user  only:  negative  angles  are  the  result  of  errors  In  the 
construction  of  the  geometry. 

3.  Dihedrals  can  normally  only  assume  definable  angles.  If  atom  i 
makes  a  dihedral  with  atoms  j,  k,  and  1,  and  the  three  atoms  j, 
k,  and  1  are  in  a  straight  line,  then  the  dihedral  has  no 
definable  angle.  During  the  calculation  this  constraint  is 
checked  continuously,  and  if  atoms  j,  k,  and  1  lie  within  0.1 
Angstroms  of  a  straight  line  the  calculation  will  output  an 
error  message  and  then  stop.  Two  exceptions  to  this  constraint 
are: 


(a)  If  the  angle  is  zero  or  180  degrees,  in  which  case  the 
dihedral  Is  not  used. 

(b)  If  atoms  j,  k,  and  1  lie  in  an  exactly  straight  line 
(usually  the  result  of  a  symmetry  constraint),  as  in  acetylene, 
acetonitrile,  but-2-yne,  etc. 


If  the  exceptions  are  used,  care  must  be  taken  to  ensure  that  the 
program  does  not  violate  these  constraints  during  any  optimizations  or 
during  any  calculations  of  derivatives  -  see  also  FORCE. 

CARTESIAN  COORDINATE  DEFINITION 

Cartesian  coordinates  consist  of  the  chemical  symbol  or  atomic 
number,  then  the  cartesian  coordinates  and  optimization  flags  but  no 
connectivity. 

MOPAC  uses  the  lack  of  connectivity  to  indicate  that  cartesian 
coordinates  are  to  be  used.  A  unique  case  is  the  triatomics  for  which 
only  Internal  coordinates  are  allowed.  This  is  to  avoid  conflict  of 
definitions:  the  user  does  not  need  to  define  the  connectivity  of  atom 
2,  and  can  elect  to  use  the  default  connectivity  for  atom  3.  As  a 
result,  a  trlatomlc  may  have  no  explicit  connectivity  defined,  the  user 
thus  taking  advantage  of  the  default  connectivity.  Since  internal 
coordinates  are  more  commonly  used  than  cartesian,  the  above  choice  was 
made. 


If  the  key-word  XYZ  is  absent  then  every  coordinate  must  be  marked 
for  optimization.  If  any  coordinates  are  not  to  be  optimized,  then  the 
key-word  XYZ  must  be  present.  The  coordinates  of  all  atoms,  Including 
atoms  1,  2  and  3  can  be  optimized.  Dummy  atoms  should  not  be  used,  for 
obvious  reasons. 
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3.2  DEFINITION  OF  ELEMENTS  AND  ISOTOPES 

Elements  are  defined  In  terms  of  their  chemical  symbols. 
Acceptable  symbols  for  MNDO  are: 

H  LI  Be  B  C  N  0  F  A1  SI  P  S  Cl  Br  I  Hg  Pb  +  -  Xx 
LI  BE  AL  SI  CL  BR  HG  PB  XX 

and  for  AMI,  acceptable  symbols  are 

H  C  N  0  F  Cl  Br  I  Xx  +  - 
CL  BR  XX 


H 

B 

C 

N 

0 

F 

SI 

P 

S 

Cl 


D1 atomics  Parameterized  within  the  MINDO/3  Formalism 


H  B 


N 


F  SI 


S  Cl 


A  star  (*)  Indicates 
that  the  atom-pair  Is 
parameterized  within 
MINDO/3. 


Extra  entitles  available  to  MNDO,  MINDO/3  and  AMI 

+  A  1001  Ionic  alkali  metal. 

-  A  1001  Ionic  halogen-like  atom 

Element  XX  Is  known  as  a  dummy  atom,  and  Is  used  In  the  definition 
of  the  geometry;  It  Is  deleted  automatically  from  any  cartesian 
coordinate  geometry  files.  Dummy  atoms  are  pure  mathematic  points,  and 
are  useful  In  defining  geometries;  for  example.  In  ammonia  the 
definition  of  C3v  symmetry  Is  facilitated  by  using  one  dummy  atom  and 
symmetry  relating  the  three  hydrogens  to  It. 

Isotopes  are  used  In  conjunction  with  chemical  symbols.  If  no 
Isotope  Is  specified,  the  average  Isotopic  mass  Is  used,  thus  chlorine 
Is  35.453. 

The  sparkles  +  and  -  have  no  mass;  If  they  are  to  be  used  In  a 
force  calculation,  then  appropriate  masses  should  be  used. 

Each  Internal  coordinate  Is  followed  by  an  Integer,  to  Indicate  the 
action  to  be  taken. 

Integer  Action 


1 

0 

-1 


Optimize  the  Internal  coordinate. 

Do  not  optimize  the  Internal  coordinate. 
Reaction  coordinate,  or  grid  Index. 
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Remarks : 

Only  one  reaction  coordinate  is  allowed,  but  this  can  be  made  more 
versatile  by  the  use  of  SYMMETRY.  If  a  reaction  coordinate  is  used,  the 
values  of  the  reaction  coordinate  should  follow  Immediately  after  the 
geometry  and  any  symmetry  data.  No  terminator  is  required,  and 
free-format-type  Input  is  acceptable. 

Using  Internal  coordinates,  the  first  atom  has  three  unoptimlzable 
coordinates,  the  second  atom  two,  (the  bond-length  can  be  optimized)  and 
the  third  atom  has  one  unoptimizable  coordinate.  None  of  these  six 
unoptimlzable  coordinates  at  the  start  of  the  geometry  should  be  marked 
for  optimization.  If  any  are  so  marked,  a  warning  is  given,  and  the 
calculation  will  continue. 

In  cartesian  coordinates  all  parameters  can  be  optimized. 


3.3  EXAMPLES  OF  COORDINATE  DEFINITIONS. 

First,  formaldehyde.  By  definition  atom  1,  here  oxygen,  is  at  the 
origin.  Atom  2  Is  defined  as  being  along  the  "x"  axis,  and  is  bonded  to 
atom  1.  Because  of  these  definitions,  the  user  does  not  need  to  specify 
the  connectivity  of  atom  2,  but  can  do  so  if  desired.  Further,  by 
default,  but  not  by  definition,  atom  3  is  connected  to  atom  2,  and  makes 
an  angle  with  atom  1.  If  the  user  wishes,  either  the  default 

connectivity  for  atom  3,  or  the  alternative  connectivity,  to  atom  1, 

making  an  angle  with  atom  2,  can  be  explicitly  defined. 

In  this  example,  the  minimum  data  necessary  to  define  the  geometry 

has  been  entered.  Thus  for  atoms  2  and  3  the  connectivity  has  been 

omitted:  the  default  connectivity  is  used.  For  atom  4,  however,  the 
full  connectivity  must  be  explicitly  defined;  no  defaults  are  allowed. 

0  The  first  atom  has  no  coordinates. 

C  1.2  1  The  C-0  bond  length  is  to  be  optimized. 

H2.0140  1.0  1.120  1  The  third  atom  is  a  Deuterium. 

H  0.0  0  0.0  1  180  0213  The  fourth  atom  is  a  Hydrogen,  and 

SYMMETRY  must  have  been  specified,  in 
order  to  not  have  to  give  a  bond-length  or  an  angle.  The  dihedral  is  point- 
group  defined  as  180  degrees. 


M0PAC  can  generate  data-files,  both  in  the  Archive  files,  and  at 
the  end  of  the  normal  output  file,  when  a  job  ends  prematurely  due  to 
time  restrictions.  Here  the  coordinate  definition  for  formic  acid  is 
shown. 
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Note  that  all  coordinates  are  generated,  as  Is  the  full 
connectivity.  Also,  the  data  are  all  neatly  lined  up.  This  Is,  of 
course,  characteristic  of  machine-generated  data,  but  Is  useful  when 
checking  for  errors. 


Format  of  internal  coordinates  in  ARCHIVE  file 


0  0.000000  0 

C  1.209615  1 

0  1.313679  1 

H  0.964468  1 

H  1.108040  1 

0  0.000000  0 


0.000000 

0 

0 

0.000000 

0 

0 

116.886168 

1 

0 

115.553316 

1 

0 

128.726078 

1 

180 

0.000000 

0 

0 

000000  0000 
000000  0100 
000000  0210 
000000  0321 
000000  0213 
000000  0000 


CHAPTER  4 


EXAMPLES 


In  this  chapter  various  examples  of  data-flles  are  described.  With 
MOPAC  comes  two  sets  of  data  for  r. inning  calculations.  One  of  these  Is 
called  MNRSD1.DAT,  and  this  will  now  be  described. 


4.1  MNRSD1  TEST  DATA  FILE  FOR  FORMALDEHYDE 

The  following  file  Is  suitable  for  generating  the  results  described 
In  the  next  section,  and  would  be  suitable  for  debugging  data. 

Line  1:  SYMMETRY 

Line  2:  Formaldehyde,  for  Demonstration  Purposes 
Line  3: 

Line  4:  0 

Line  5:  C  1.2  1 

Line  6:  H  1.1  1  120  1 

Line  7:  H  1 . 1  0  120  0  180  0  2  1  3 

Line  8: 

Line  9:  314 

Line  10:  324 

Line  11: 


This  data  could  be  more  neatly  written  as 


Ltne 

Line 

Line 

1 

2 

3 

SYMMETRY 

Formaldehyde,  for 

Demonstration  Purposes 

Line 

4 

0 

0.000000 

0 

0.000000 

0 

0.000000 

0 

0 

0 

0 

Line 

5 

C 

1 .200000 

1 

0.000000 

0 

0.000000 

0 

1 

0 

0 

Line 

6 

H 

1.100000 

1 

120.000000 

1 

0.000000 

0 

2 

1 

0 

Line 

7 

H 

1.100000 

0 

120.000000 

0 

180.000000 

0 

2 

1 

3 

Line 

Line 

8 

9 

0 

3. 

0.000000 

1.  A, 

0 

0.000000 

0 

0.000000 

0 

0 

0 

0 
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These  two  data-flles  will  produce  Identical  results  files. 

In  all  geometric  specifications  care  must  be  taken  In  defining  the 
Internal  coordinates  to  ensure  that  no  three  atoms  being  used  to  define 
a  fourth  atom's  dihedral  angle  ever  fall  Into  a  straight  line.  This  can 
happen  In  the  course  of  a  geometry  optimization.  In  a  SADDLE  calculation 
or  In  following  a  reaction  coordinate.  If  such  a  condition  should 
develop,  then  the  position  of  the  dependent  atom  would  become 
Ill-defined. 


r. 


4.2  MOPAC  OUTPUT  FOR  TEST-DATA  FILE  MNRSD1 


***  FRANK  J  SEILER  RES.  LAB.,  U.S.  AIR  FORCE  ACADEMY,  COLO.  SPGS.,  CO.  80840  *** 

MNDO  CALCULATION  RESULTS  Note  1 


*  VERSION  3.10-PC  Note  2 

*  SYMMETRY  -  SYMMETRY  CONDITIONS  TO  BE  IMPOSED 


PARAMETER  DEPENDENCE  DATA 

REFERENCE  ATOM  FUNCTION  NO.  DEPENDENT  ATOM(S) 

3  1  4 

3  2  4 

DESCRIPTIONS  OF  THE  FUNCTIONS  USED 

1  BOND  LENGTH  IS  SET  EQUAL  TO  THE  REFERENCE  BOND  LENGTH 

2  BOND  ANGLE  IS  SET  EQUAL  TO  THE  REFERENCE  BOND  ANGLE 

SYMMETRY  Note  3 

Formaldehyde,  for  Demonstration  Purposes 


ATOM 

CHEMICAL 

BOND  LENGTH 

BOND  ANGLE 

TWIST  ANGLE 

NUMBER 

SYMBOL 

(ANGSTROMS) 

(DEGREES) 

(DEGREES) 

(I) 

NA:  I 

NB:NA:I 

NC:NB:NA:I 

NA 

NB 

NC 

1 

0 

Note 

4 

2 

C 

1 .20000  * 

1 

3 

H 

1.10000  * 

120.00000  * 

2 

1 

4 

H 

1.10000 

120.00000 

180.00000 

2 

1 

3 

CARTESIAN  COORDINATES 

NO.  ATOM  X  Y  Z 
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1 

2 

3 

4 

H:  (MNDO) 
C:  (MNDO) 
0:  (MNDO) 


0  0.0000  0.0000  0.0000 

C  1.2000  0.0000  0.0000  Note  5 

H  1.7500  0.9526  0.0000 

H  1.7500  -0.9526  0.0000 

M.J.S.  DEWAR,  W.  THIEL,  0.  AM.  CHEM.  SOC. ,  99,  4907,  (1977) 

M.J.S.  DEWAR,  W.  THIEL,  0.  AM.  CHEM.  SOC.,  99,  4907,  (1977) 

M.J.S.  DEWAR,  W.  THIEL,  J.  AM.  CHEM.  SOC.,  99,  4907,  (1977) 


RHF  CALCULATION,  NO.  OF  DOUBLY  OCCUPIED  LEVELS  =  6 


INTERATOMIC  DISTANCES 

0  1  C  2  H  3  H  4 


0  1  0.000000 

C  2  1.200000  0.000000 

H  3  1.992486  1.100000  0.000000 

H  4  1.992486  1.100000  1.905256  0. 000000 

CYCLE:  1  TIME:  43.00  TIME  LEFT:  3522.0  GRAD.: 

CYCLE:  2  TIME:  17.00  TIME  LEFT:  3505.0  GRAD.: 

HEAT  OF  FORMATION  TEST  SATISFIED 
PETERS  TEST  SATISFIED 
SYMMETRY 

Formaldehyde,  for  Demonstration  Purposes 


Note  6 


36.370  HEAT:  -32.86455 
8.250  HEAT:  -32.88052 

Note  7 
Note  8 
Note  9 
Note  10 


PETERS  TEST  WAS  SATISFIED  IN  FLETCHER-POWELL  OPTIMISATION  Note  11 

SCF  FIELD  WAS  ACHIEVED  Note  12 


MNDO 

CALCULATION 

Note  13 

VERSION  3.10-PC 

FINAL  HEAT  OF  FORMATION 

3 

-32.88189 

KCAL 

Note  14 

ELECTRONIC  ENERGY 

-870.74066 

EV 

CORE-CORE  REPULSION 

3 

392.62149 

EV 

IONISATION  POTENTIAL 

= 

11.04150 

06-03-86 

NO.  OF  FILLED  LEVELS 

= 

6 

MOLECULAR  WEIGHT 

= 

30.026 

SCF  CALCULATIONS 
COMPUTATION  TIME 


15 

127.00'  SECONDS 


Note  15 
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ATOM 

CHEMICAL 

BONO  LENGTH 

BOND  ANGLE 

TWIST  ANGLE 

NUMBER 

SYMBOL 

(ANGSTROMS) 

(DEGREES) 

(DEGREES) 

(I) 

NA:  I 

NB:NA:I 

NC : NB : NA : I 

NA  NB  NC 

1 

0 

Note  16 

2 

C 

1 .21646  * 

1 

3 

H 

1.10601  * 

123.51050  * 

2  1 

4 

H 

1.10601 

123.51050 

180.00000 

2  1  3 

INTERATOMIC 

DISTANCES 

0  1 

C  2 

H  3 

H  4 

0  1  0.000000 

C  2  1.216455  0.000000 

H  3  2.046609  1.106012  0.000000 

H  4  2.046609  1.106012  1.844352  0.000000 


EIGENVALUES 


-42.99227  -25.11986 

-16.95504  - 

■16.29961  - 

14.17985  -11.04150 

0.85990 

3.67652 

3.84972 

7.12622 

Note 

17 

NET  ATOMIC  CHARGES  AND  DIPOLE  CONTRIBUTIONS 

ATOM 

1  NO.  TYPE 

CHARGE 

ATOM  ELECTRON 

DENSITY 

1 

0 

-0.2902 

6.2902 

2 

C 

0.2922 

3.7078 

Note 

18 

3 

H 

-0.0010 

1.0010 

4 

H 

-0.0010 

1.0010 

DIPOLE 

X 

Y 

Z 

TOTAL 

POINT-CHG. 

1.690 

0.000 

0.000 

1.690 

HYBRID 

0.475 

0.000 

0.000 

0.475 

Note 

19 

SUM 

2.164 

0.000 

0.000 

2.164 

CARTESIAN  COORDINATES 

NO. 

ATOM 

X 

Y  Z 

1 

0 

0.0000 

0.0000  0.0000 

2 

C 

1 .2165 

0.0000  0.0000 

3 

H 

1 .8271 

0.9222  0.0000 

4 

H 

1 .8271 

-0.9222  0.0000 

ATOMIC  ORBITAL  ELECTRON  POPULATIONS 

1.88260 

1.21603 

1 .89109 

1.30047 

1.25526  0.86214 

0.89092 

0.69953 

1.00098 

1.00098 

Note 

20 
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NOTES  ON  RESULTS  FILE 

NOTE  1:  The  banner  indicates  whether  the  calculation  uses  a  MNDO, 
MINDO/3  or  AMI  Hamiltonian;  here,  the  default  MNDO  Hamiltonian  Is  used. 

NOTE  2:  The  Version  number  Is  a  constant  for  any  release  of  MOPAC, 
and  refers  to  the  program,  not  to  the  Hamiltonians  used.  The  version 
number  should  be  cited  In  any  correspondence  regarding  MOPAC. 

All  the  key-words  used,  along  with  a  brief  explanation,  should  be 
printed  at  this  time.  If  a  key-word  Is  not  printed.  It  has  not  been 
recognized  by  the  program.  Key-words  can  be  in  upper  or  lower  case 
letters,  or  any  mixture.  Originally  only  uppercase  letters  were 
allowed,  as  use  of  lower  case  letters  involves  an  extension  to  the 
FORTRAN-77  standard.  It  was  found,  however,  that  all  computers  used 
support  lower  case  letters,  so  the  ability  to  use  lower  case  letters  has 
been  included,  for  the  sake  of  convenience. 

NOTE  3:  Symmetry  information  Is  output  to  allow  the  user  to  verify 
that  the  requested  symmetry  functions  have  In  fact  been  recognized  and 
used. 

NOTE  4:  The  data  for  this  example  used  a  mixture  of  atomic  numbers 
and  chemical  symbols,  but  the  Internal  coordinate  output  is  consistently 
in  chemical  symbols. 

The  atoms  In  the  system  are,  In  order: 

Atom  1,  an  oxygen  atom;  this  is  defined  as  being  at  the 
origin. 

Atom  2,  the  carbon  atom.  Defined  as  being  1.2  Angstroms 
from  the  oxygen  atom.  It  Is  located  In  the  +x  direction.  This 
distance  is  marked  for  optimization. 

Atom  3,  a  hydrogen  atom.  It  Is  defined  as  being  1.1 
Angstroms  from  the  carbon  atom,  and  making  an  angle  of  120 
degrees,  with  the  oxygen  atom.  The  asterisks  Indicate  that  the 
bond  length  and  angle  are  both  to  be  optimized. 

Atom  4,  a  hydrogen  atom.  The  bond  length  supplied  has  been 
overwritten  with  the  symmetry-defined  C-H  bond  length.  Atom  4  Is 
defined  as  being  1.1  Angstroms  from  atom  2,  making  a  bond-angle 
of  120  degrees  with  atom  1,  and  a  dihedral  angle  of  180  degrees 
with  atom  3. 

None  of  the  coordinates  of  atom  4  are  marked  for 
optimization.  The  bond-length  and  angle  are  symmetry-defined  by 
atom  3,  and  the  dihedral  Is  group-theory  symmetry-defined  as 
being  180  degrees.  (The  molecule  Is  flat.) 

NOTE  5:  The  cartesian  coordinates  are  calculated  as  follows: 
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Stage  1:  The  coordinate  of  the  first  atom  Is  defined  as 
being  at  the  origin  of  cartesian  space,  while  the  coordinate  of 
the  second  atom  Is  defined  as  being  displaced  by  Its  defined  bond 
length  along  the  positive  x-axIs.  The  coordinate  of  the  third 
atom  Is  defined  as  being  displaced  by  Its  bond  length  In  the  x-y 
plane,  from  either  atom  1  or  2  as  defined  In  the  data,  or  from 
atom  2  If  no  numbering  Is  given.  The  angle  It  makes  with  atoms  1 
and  2  Is  that  given  by  Its  bond  angle. 

The  dihedral,  which  first  appears  In  the  fourth  atom.  Is 
defined  according  to  the  I.U.P.A.C.  convention.  NOTE:  This  Is 
different  from  previous  versions  of  MNDO  and  MINDO/3,  where  the 
dihedral  had  the  opposite  chirality  to  that  defined  by  the 
I.U.P.A.C.  convention. 

Stage  2:  Any  dummy  atoms  are  removed.  As  this  particular 
system  contains  no  dummy  atoms,  nothing  Is  done. 

NOTE  6:  The  Interatomic  distances  are  output  for  the  user's 
advice,  and  a  simple  check  made  to  Insure  that  the  smallest  Interatomic 
distance  Is  greater  than  0.8  Angstroms. 

NOTE  7:  The  geometry  Is  optimized  In  a  series  of  cycles,  each 
cycle  consisting  of  a  line-search  and  calculation  of  the  gradients.  The 
time  given  Is  the  C.P.U.  time  for  the  cycle;  time  left  Is  the  total 
time  requested  (here  100  seconds)  less  the  C.P.U.  time  since  the  start 
of  the  calculation  (which  Is  earlier  than  the  start  of  the  first 
cycle!).  These  times  can  vary  considerably  from  cycle  to  cycle  due  to 
different  options  being  used,  for  example  whether  or  not  the  density 
matrix  Is  being  updated  at  every  point  In  the  line  search.  The  gradient 
Is  the  scalar  length  In  kcal /mole /Angstrom  of  the  gradient  vector. 

NOTE  8:  At  the  end  of  the  Davldon-Fletcher-Powel 1  geometry 
optimization  a  message  Is  given  which  Indicates  how  the  optimization 
ended.  All  "normal"  termination  messages  contain  the  word  "satisfied"; 
other  terminations  may  give  acceptable  results,  but  more  care  should  be 
taken,  particularly  regarding  the  gradient  vector. 

NOTE  9  and  10:  The  key-words  used,  titles  and  comments  are 
reproduced  here  to  remind  the  user  of  the  name  of  the  calculation. 

NOTE  11  and  12:  Two  messages  are  given  here.  The  first  Is  a 
reminder  of  how  the  geometry  was  obtained,  whether  from  the  D.F.P.  or 
Bartel's  method.  For  any  further  results  to  be  printed  the  second 
message  must  be  as  shown;  when  no  SCF  Is  obtained  no  results  will  be 
printed. 

NOTE  13:  Again,  the  results  are  headed  with  either  MNDO  or  MINDO/3 
banners,  and  the  version  number. 

NOTE  14:  The  total  energy  of  the  system  Is  the  addition  of  the 
electronic  and  nuclear  terms.  The  heat  of  formation  Is  relative  to  the 
elements  In  their  standard  state.  The  I.P.  Is  the  negative  of  the 
energy  level  of  the  highest  occupied,  or  highest  partial ly  occupied 
molecular  orbital  (In  accordance  with  Koopmans'  theorem). 
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NOTE  15:  Advice  on  time  required  for  the  calculation.  This  Is 
obviously  useful  In  estimating  the  times  required  for  other  systems. 

NOTE  16:  The  fully  optimized  geometry  Is  printed  here.  If  a 
parameter  Is  not  marked  for  optimization.  It  will  not  be  changed  unless 
It  Is  a  symmetry-related  parameter. 

NOTE  17:  The  roots  are  the  eigenvalues  or  energy  levels  In 
electron  volts  of  the  molecular  orbitals.  There  are  six  filled  levels, 
therefore  the  HOMO  has  an  energy  of  -11. 04? eV,  and  Is  mainly  lone-pair 
on  oxygen.  The  eigenvectors  form  an  orthonormal  set. 

NOTE  18:  The  charge  on  an  atom  is  the  sum  of  the  positive  core 
charge;  for  hydrogen,  carbon,  and  oxygen  these  numbers  are  1.0,  4.0,  and 
6.0,  respectively,  and  the  negative  of  the  number  of  valence  electrons, 
or  atom  electron  density  on  the  atom,  here  1.0010,  3.7079,  and  6.2902 
respectively. 

NOTE  19:  The  dipole  Is  the  scalar  of  the  dipole  vector  in 
cartesian  coordinates.  The  components  of  the  vector  coefficients  are 
the  point-charge  dipole  and  the  hybridization  dipole.  In  formaldehyde 
there  Is  no  z-dlpole  since  the  molecule  Is  flat. 

NOTE  20:  MND0  AMI,  and  MIND0/3  all  use  the  Coulson  density  matrix. 
Only  the  diagonal  elements  of  the  matrix,  representing  the  valence 
orbital  electron  populations,  will  be  printed,  unless  the  keyword 
DENSITY  Is  specified. 


CHAPTER  5 


TESTDATA 


This  example  Is  taken  from  the  first  data-flle  In  TESTDATA.DAT,  and 
Illustrates  the  working  of  a  FORCE  calculation. 


5.1 

DATA  FILE 

FOR  A  FORCE  CALCULATION 

Line 

1 

ROT-2 

THERMOC298 

,298,)  PRECISE 

FORCE  ISOTOPE 

SYMMETRY 

Line 

2 

DEMONSTRATION  OF  MOPAC  -  FORCE  AND  THERMODYNAMICS  CALCULATION 

Line 

3 

FORMALDEHYDE, 

MNDO  ENERGY  -  - 

-32. 

8819 

Line 

4a 

0 

0.000000 

0  0.000000 

0 

0.000000 

0 

0 

0 

0 

Line 

4b 

C 

1.216486 

1  0.000000 

0 

0.000000 

0 

1 

0 

0 

Line 

4c 

H 

1.106188 

1  123.515021 

1 

0.000000 

0 

2 

1 

0 

Line 

4d 

H 

1.106188 

1  123.515021 

1 

179.995690 

1 

2 

1 

3 

Line 

4e 

0 

0.000000 

0  0.000000 

0 

0.000000 

0 

0 

0 

0 

Line 

5a 

3, 

1,  4, 

Line 

5b 

3, 

2,  4, 

5.2  RESULTS  FILE  FOR  THE  FORCE  CALCULATION 


***  FRANK  J  SEILER  RES.  LAB.,  U.S.  AIR  FORCE  ACADEMY,  COLO.  SPGS.,  CO.  80840  *** 


MNDO  CALCULATION  RESULTS 


* 


* 


VERSION  3.10-PC 

SYMMETRY  -  SYMMETRY  CONDITIONS  TO  BE  IMPOSED 
FORCE  -  FORCE  CALCULATION  SPECIFIED 
PRECISE  -  CRITERIA  TO  BE  INCREASED  BY  100  TIMES 
ISOTOPE  -  FORCE  MATRIX  WRITTEN  TO  DISK  (CHAN.  9  ) 
THERMO  -  THERMODYNAMIC  QUANTITIES  TO  BE  CALCULATED 
ROT  '  -  SYMMETRY  NUMBER  OF  2  SPECIFIED 
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PARAMETER  DEPENDENCE  DATA 

REFERENCE  ATOM  FUNCTION  NO.  DEPENDENT  ATOM(S) 
3  1  4 

3  2  4 


DESCRIPTIONS  OF  THE  FUNCTIONS  USED 

1  BOND  LENGTH  IS  SET  EQUAL  TO  THE  REFERENCE  BOND  LENGTH 

2  BOND  ANGLE  IS  SET  EQUAL  TO  THE  REFERENCE  BOND  ANGLE 

ROT=2  THERMO(298,298, )  PRECISE  FORCE  ISOTOPE  SYMMETRY  Note  1 

DEMONSTRATION  OF  MOPAC  -  FORCE  AND  THERMODYNAMICS  CALCULATION 
FORMALDEHYDE,  MNDO  ENERGY  *  -32.8819 


ATOM 

CHEMICAL 

BOND  LENGTH 

BOND  ANGLE 

TWIST  ANGLE 

NUMBER 

SYMBOL 

(ANGSTROMS) 

(DEGREES) 

(DEGREES) 

(I) 

NA:  I 

NB:NA:I 

NC : NB : NA : I 

NA 

NB 

1 

0 

2 

C 

1 .21649  * 

1 

3 

H 

1.10619  * 

123.51502  * 

2 

1 

4 

H 

1.10619  * 

123.51502  * 

179.99569  * 

2 

1 

CARTESIAN 

COORDINATES 

- 

NO. 

ATOM 

X 

Y 

Z 

1 

0 

0.0000 

0.0000 

0.0000 

2 

C 

1.2165 

0.0000 

0.0000 

3 

H 

1 .8273 

0.9223 

0.0000 

4 

H 

1.8273 

-0.9223 

0.0001 

(MNDO):  M.J.S. 

,  DEWAR,  W. 

THIEL,  J. 

AM.  CHEM.  SOC.,  99,  4907,  (1977) 

C:  (MNDO):  M.J.S.  DEWAR,  W.  THIEL,  J.  AM.  CHEM.  SOC. ,  99,  4907,  (1977) 

0:  (MNDO):  M.J.S.  DEWAR,  W.  THIEL,  J.  AM.  CHEM.  SOC.,  99,  4907,  (1977) 

RHF  CALCULATION,  NO.  OF  DOUBLY  OCCUPIED  LEVELS  =  6 


INTERATOMIC  DISTANCES 

0  1  C  2  H  3  H  4 


0  1  0.000000 

C  2  1.216486  0.000000 

H  3  2.046832  1.106188  0.000000 

H  4  2.046832  1.106188  1.844549  0.000000 


HEAT  OF  FORMATION  =  -32.881894  KCALS/MOLE 


-  44  - 


TESIUATA 
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INTERNAL  COORDINATE  DERIVATIVES 


ATOM  AT.  NO.  BOND 


ANGLE  DIHEDRAL 


1  0 

2  C  0.020543 

3  H  0.140798  0.010501 

4  H  0.000000  0.000000  -0.002392 


GRADIENT  NORM  -  0.14270 


TIME  FOR  SCF  CALCULATION  -  17.0 

TIME  FOR  DERIVATIVES  -  18.0 


SYMMETRY  WAS  SPECIFIED,  BUT  CANNOT  BE  USED  HERE 
MOLECULAR  WEIGHT  -  30.03 


Note  2 


Note  3 


PRINCIPAL  MOMENTS  OF  INERTIA  IN  CM(-l) 

A-  9.831716  B  =  1.261967  C-  1.118411 


PRINCIPAL  MOMENTS  OF  INERTIA  IN  UNITS  OF  10**(-40)*GRAM-CM*‘2 


A-  2.847177  B  -  22.181748  C-  25.028925 


ORIENTATION  OF  MOLECULE  IN  FORCE  CALCULATION 


NO.  ATOM 


Z 


1 

2 

3 

4 


8  -0.6093  0.0000 
6  0.6072  0.0000 
1  1.2180  0.9223 
1  1.2180  -0.9223 


0.0000 

0.0000 

0.0000 

0.0001 


FIRST  DERIVATIVES  WILL  BE  USED  IN  THE  CALCULATION  OF  SECOND  DERIVATIVES 
DEFAULT  TIME  OF  3600.00  SECONDS  ALLOCATED  FOR  THIS  STEP 


ESTIMATED  TIME  TO  COMPLETE  CALCULATION  -  840.00  SECONDS 


STEP: 

1  TIME  - 

46.00 

SECS, 

INTEGRAL  - 

46.00 

TIME 

LEFT: 

3519.00 

STEP: 

2  TIME  - 

50.00 

SECS, 

INTEGRAL  - 

96.00 

TIME 

LEFT: 

3469.00 

STEP: 

3  TIME  - 

49.00 

SECS, 

INTEGRAL  - 

145.00 

TIME 

LEFT: 

3420.00 

STEP: 

4  TIME  - 

47.00 

SECS, 

INTEGRAL  - 

192.00 

TIME 

LEFT: 

3373.00 

sTEP: 

5  TIME  - 

49.00 

SECS, 

INTEGRAL  - 

241 .00 

TIME 

LEFT: 

3324.00 
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STEP 

6 

TIME 

— 

47.00 

SECS, 

INTEGRAL 

= 

288.00 

TIME 

LEFT 

3277.00 

STEP 

7 

TIME 

3 

46.00 

SECS, 

INTEGRAL 

a 

334.00 

TIME 

LEFT 

3231 .00 

STEP 

8 

TIME 

3 

45.00 

SECS, 

INTEGRAL 

= 

379.00 

TIME 

LEFT 

3186.00 

STEP 

9 

TIME 

3 

53.00 

SECS, 

INTEGRAL 

a 

432.00 

TIME 

LEFT 

3133.00 

STEP 

10 

TIME 

, 

48.00 

SECS, 

INTEGRAL 

= 

480.00 

TIME 

LEFT 

3085.00 

STEP 

11 

TIME 

3 

45.00 

SECS, 

INTEGRAL 

a 

525.00 

TIME 

LEFT 

3040.00 

STEP 

12 

TIME 

3 

52.00 

SECS, 

INTEGRAL 

a 

577.00 

TIME 

LEFT 

2988.00 

FORCE  MATRIX  IN  MILLIDYNES/ANGSTROM 

0  1  C  2  H  3  H  4 


0  1 
C  2 
H  3 
H  4 


9.557769 

8.683380  11.426325 

0.598735  2.552483  3.033980 

0.598769  2.552451  0.303925  3.033952 


HEAT  OF  FORMATION  =  -32.881894  KCALS/MOLE 


TRIVIAL  VIBRATIONS,  SHOULD  BE  ZERO 
0 . 0000=TX  0.0000=TY  0.0000=TZ  0.0022-RX 


Note  4 

0.0003=RY  0.0001=RZ 


FORCE  CONSTANTS  IN  MILLIDYNES/ANGSTROM  (=  10**5  DYNES/CM) 
0.94406  1.35190  1.86920  4.07261  6.74514  18.72184 

ASSOCIATED  EIGENVECTORS 


Note  5 


NO, 

1  ■ 

2 

3 

4 

5 

6 

0.94406 

1.35190 

1.86920 

4.07261 

6.74514 

18.72184 

1 

-0.33016 

0.00000 

-0.00002 

0.41812 

0.00000 

-0.68277 

2 

0.00000 

0.60834 

-0.00002 

0.00000 

-0.04811 

0.00000 

3 

-0.00001 

-0.00003 

-0.28948 

0.00000 

0.00000 

0.00000 

4 

-0.34275 

0.00000 

-0.00002 

0.32022 

0.00000 

0.72800 

5 

0.00000 

-0.34892 

0 . 00004 

0.00000 

0.79262 

0.00000 

6 

-0.00001 

0.00003 

0.86603 

0.00004 

-0.00003 

0.00000 

7 

0.33646 

-0.48710 

0.00003 

-0.36917 

-0.21481 

-0.02261 

8 

-0.52302 

-0.12971 

0.00001 

-0.47438 

-0.37226 

-0.03759 

9 

0.00003 

-0.00001 

-0.28827 

0.00000 

0.00000 

0.00000 

10 

0.33645 

■0.48710 

0.00001 

-0.36917 

0.21480 

-0.02261 

11 

0.52303 

-0.12971 

-0.00003 

0.47438 

-0.37226 

0.03759 

12 

-0.00001 

0.00000 

-0.28827 

-0.00004 

0.00003 

0.00000 
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I 


ZERO  POINT  ENERGY  18.150  KILOCALORIES  PER  MOLE 


THE  LAST  6  VIBRATIONS  ARE  THE  TRANSLATION  AND  ROTATION  MODES 

THE  FIRST  THREE  OF  THESE  BEING  TRANSLATIONS  IN  X,  Y,  AND  Z,  RESPECTIVELY 


FREQUENCIES,  REDUCED  MASSES  AND  VIBRATIONAL  DIPOLES 


R 

i 

1 

2 

3 

4 

5 

6 

S 

FREQ(I) 

1209.9907 

1215.0493 

1490.5916 

2114.5759 

3255.3657 

3301.5822 

MASS(I) 

1 .93736 

1 . 50020 

0.66394 

6.79206 

0.75075 

0.66443 

■ 

DIPX(I) 

-0.00002 

0.00033 

-0.34444 

-3.36646 

-0.00009 

-0.34788 

DIPY(I) 

-0.85473 

-0.00087 

-0.00006 

0.00002 

0.78304 

-0.00016 

ii 

DIPZ(I) 

-0.00010 

0.12739 

-0.00002 

-0.00008 

-0.00003 

-0.00001 

t 

DIPT  < I ) 

0.85473 

0.12739 

0.34444 

3.36646 

0.78304 

0.34788 

$ 

I 

7 

8 

9 

10 

11 

12 

V 

FREQ(I) 

-0.0056 

0.0005 

-0.0019 

87.0251 

11 .1816 

10.6530 

< 

*„* 

MASS(I) 

0.00000 

0.00000 

0.00000 

0.50395 

0.00000 

0.00000 

DIPX(I) 

0.00002 

0.00001 

0.00023 

0.00000 

-0.00028 

-0.00001 

■ 

DIPY(I) 

-0.00004 

-0.00005 

-0.00005 

0.00000 

0.00005 

0.46525 

[3 

DIPZ(I) 

0.00000 

0.00000 

0.00004 

0.00001 

0.56207 

-0.00002 

k: 

DIPT(I) 

0.00005 

0.00005 

0.00024 

0.00002 

0.56207 

0.46525 

Note  6 


NORMAL  VECTORS 


Note  7 


k 

ROOT  NO. 

1 

2 

3 

4 

5 

6 

■ 

1209.99070 

1215.04928 

1490.59161 

2114.57589 

3255.36573 

3301 .58223 

s 

1 

0.00000 

0.00000 

-0.04332 

0.52497 

0.00000 

0.00070 

2 

0.07922 

0.00007 

0.00000 

0.00000 

0.00432 

0.00000 

3 

0.00004 

-0.04399 

0.00000 

0.00000 

0.00000 

0.00000 

U 

4 

0.00000 

0.00000 

-0.03977 

-0.66818 

-0.00001 

-0.06524 

5 

-0.15855 

-0.00015 

0.00000 

0 . 00000 

0.08917 

-0.00002 

*  * 

6 

-0.00016 

0.17532 

0.00000 

0.00001 

0.00000 

0.00000 

r, 

* 

7 

-0.62003 

-0.00060 

0.58074 

-0.18535 

-0.41965 

0.38325 

C | 

8 

0.31599 

0.00030 

-0.40126 

-0.32346 

-0.56544 

0.59260 

£ 

9 

0.00067 

-0.69546 

0.00001 

-0.00002 

0.00000 

-0.00001 

10 

0.62003 

0.00060 

0.58074 

-0.18533 

0.41979 

0.38310 

F 

11 

0.31600 

0.00025 

0.40126 

0.32344 

-0.56567 

-0.59239 

*"«£  '  -r  <  -  *. 
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12 

0.00065 

-0.69546 

-0.00002 

-0.00004 

0.00005 

0.00004 

ROOT  NO. 

7 

8 

9 

10 

11 

12 

-0.00557 

0.00049 

-0.00194 

87.02506 

11.18157 

10.65295 

1 

0.50000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

2 

0.00000 

0.50000 

0.00000 

0.00000 

-0.00001 

-0.26200 

3 

0.00000 

0.00000 

-0.50000 

0.00000 

-0.31646 

0.00001 

4 

0.50000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

5 

0.00000 

0.50000 

0.00000 

0.00000 

0.00001 

0.26110 

6 

0.00000 

0.00000 

-0.50000 

0.00000 

0.31537 

-0.00001 

7 

0.50000 

0.00000 

0.00000 

0.00000 

-0.00002 

-0.39659 

8 

0.00000 

0.50000 

0.00000 

0.00001 

0.00002 

0.52375 

9 

0.00000 

0.00000 

-0.50000 

0.70711 

0.63261 

-0.00002 

10 

0.50000 

0.00000 

0 . 00000 

0.00000 

-0.00002 

0.39659 

11 

0.00000 

0.50000 

0.00000 

-0.00005 

0.00002 

0.52375 

12 

0.00000 

0.00000 

-0.50000 

-0.70711 

0.63261 

-0.00002 

DESCRIPTION  OF  VIBRATIONS 


VIB.  1 

ATOMS 

C 

2 

AND 

H 

3 

SHIFT 

0.78 

ANGSTROMS 

6.8% 

RADIALLY 

FREQ. 

1209.99 

C 

2 

H 

4 

0.78 

6.8% 

VIB.  2 

ATOMS 

C 

2 

AND 

H 

3 

SHIFT 

0.87 

ANGSTROMS 

0.0% 

RADIALLY 

FREQ. 

1215.05 

C 

2 

H 

4 

0.87 

0.0% 

VIB.  3 

ATOMS 

C 

2 

AND 

H 

3 

SHIFT 

0.74 

ANGSTROMS 

1.1% 

RADIALLY 

FREQ. 

1490.59 

C 

2 

H 

4 

0.74 

1.1% 

VIB.  4 

ATOMS 

0 

1 

AND 

C 

2 

SHIFT 

1.19 

ANGSTROMS 

100.0% 

RADIALLY 

FREQ. 

2114.58 

C 

2 

H 

3 

0.58 

0.5% 

- 

C 

2 

H 

4 

0.58 

0.5% 

VIB.  5 

ATOMS 

c 

2 

AND 

H 

3 

SHIFT 

0.78 

ANGSTROMS 

100.0% 

RADIALLY 

FREQ. 

3255.37 

c 

2 

H 

4 

0.78 

100.0% 

VIB.  6 

ATOMS 

c 

2 

AND 

H 

3 

SHIFT 

0.74 

ANGSTROMS 

99.8% 

RADIALLY 

FREQ. 

3301.58 

c 

2 

H 

4 

0.74 

99.8% 

VIB.  10 

ATOMS 

c 

2 

AND 

H 

3 

SHIFT 

0.71 

ANGSTROMS 

0.0% 

RADIALLY 

FREQ. 

87.03 

c 

2 

H 

4 

0.71 

0.0% 

VIB.  11 

ATOMS 

0 

1 

AND 

C 

2 

SHIFT 

0.63 

ANGSTROMS 

0.0% 

RADIALLY 

FREQ. 

11.18 

c 

2 

H 

3 

0.32 

0.0% 

c 

2 

H 

4 

0.32 

0.0% 

VIB.  12 

ATOMS 

0 

1 

AND 

C 

2 

SHIFT 

0.52 

ANGSTROMS 

0.0% 

RADIALLY 

FREQ. 

10.65 

c 

2 

H 

3 

0.48 

0.0% 

SYSTEM  IS  A  GROUND  STATE 


FORMALDEHYDE,  MNDO  ENERGY  -  -32.8819 

DEMONSTRATION  OF  MOPAC  -  FORCE  AND  THERMODYNAMICS  CALCULATION 


MOLECULE  IS  NOT  LINEAR 

THERE  ARE  6  GENUINE  VIBRATIONS  IN  THIS  SYSTEM 
THIS  THERMODYNAMICS  CALCULATION  IS  LIMITED  TO 
MOLECULES  WHICH  HAVE  NO  INTERNAL  ROTATIONS 


Note  8 


CALCULATED  THERMODYNAMIC  PROPERTIES 


TEMP. 

<K) 

PARTITION  FUNCTION 

ENTHALPY 

CAL /MOL 

HEAT  CAPACITY 
CAL/K/MOL 

ENTROPY 

CAL/K/MOL 

298 

VIB. 

1 .007 

23.37536292 

0.47807894 

0.09142732 

ROT. 

709. 

888.305 

2.981 

16.026 

INT. 

714. 

911.681 

3.459 

16.117 

TRA. 

TOT. 

0. 1 59E+27 

1480.509 

2392.1893 

4.968 

8.4271 

36.113 

52.2301 

NOTE  1:  All  three  words,  ROT,  FORCE,  and  THERMO  are  necessary  In 
order  to  obtain  thermodynamic  properties.  In  order  to  obtain  results 
for  only  one  temperature,  THERMO  has  the  first  and  second  arguments 
Identical.  The  symmetry  number  for  the  C2v  point-group  Is  2. 

NOTE  2:  Internal  coordinate  derivatives  are  In  Kcal /Angstrom  or 
Kcal/radlan.  Values  of  less  than  about  0.2  are  quite  acceptable. 

NOTE  3:  In  larger  calculations,  the  time  estimates  are  useful.  In 
practice  they  are  pessimistic,  and  only  about  70%  of  the  time  estimated 
will  be  used,  usually.  The  principal  moments  of  inertia  can  be  directly 
related  to  the  microwave  spectrum  of  the  molecule.  They  are  simple 
functions  of  the  geometry  of  the  system,  and  are  usually  predicted  with 
very  high  accuracy. 

NOTE  4:  The  trivial  vibrations  are,  In  order,  translation  In  x,  y, 
and  z,  and  rotation  In  x,  y,  and  z.  Errors  In  the  force  calculation 
appear  In  Tx,  Ty,  and  Tz,  and  in  the  geometry  in  Rx,  Ry,  and  Rz. 
Ideally,  all  six  should  be  zero,  but  values  of  Tx,  Ty,  and  Tz  less  than 
about  0.0050,  and  Rx,  Ry,  and  Rz  less  than  about  0.0080  are  acceptable. 
For  linear  systems  there  will  only  be  5  trivial  vibrations,  and  for 
polymers  only  4.  For  polymers,  only  the  three  translations  are 
considered  as  trivial;  the  single  trivial  rotation  Is  considered  in  the 
genuine  force-constants,  and  should  be  very  small. 
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NOTE  5:  These  are  the  "genuine"  force  constants:  there  will 
normally  be  3N-6  of  these. 

NOTE  6:  The  transition  dipoles  are  useful  In  determining  the 
Intensity  of  I.R.  absorption  bands.  The  reduced  mass  Indicates  which 
atoms  are  Involved  In  a  vibration;  thus,  vibration  10  principally 
Involves  carbon  and  oxygen. 

NOTE  7:  These  are  the  normal  modes  of  vibration.  The  last  three, 
five  or  six  are  the  "trivial"  vibrations,  and  should  be  small. 

NOTE  8:  In  order,  the  thermodynamic  quantities  calculated  are: 
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5.3  EXAMPLE  OF  REACTION  PATH  WITH  SYMMETRY 

In  this  example,  one  methyl  group  In  ethane  Is  rotated  relative  to 
the  other  and  the  geometry  Is  optimized  at  each  point.  As  the  reaction 
coordinate  Involves  three  Hydrogen  atoms  moving,  symmetry  Is  Imposed  to 


ensure 

equivalence 

of  all  hydrogens. 

Line 

1 

SYMMETRY 

T-600 

Line 

2 

ROTATION  OF  METHYL  GROUP  IN 

ETHANE 

Line 

3 

EXAMPLE  OF  A  REACTION  PATH  CALCULATION 

Line 

4 

C 

Line 

5 

C 

1.479146  1 

Line 

6 

H 

1.109475  1 

111.328433 

1 

Line 

7 

H 

1.109470  0 

111.753160 

0 

120.000000 

0 

2 

1 

3 

Line 

8 

H 

1.109843  0 

110.103163 

0 

240.000000 

0 

2 

1 

3 

Line 

9 

H 

1.082055  0 

121 .214083 

0 

60.000000 

-1 

1 

2 

3 

Line 

10 

H 

1.081797  0 

121.521232 

0 

180.000000 

0 

1 

2 

3 

Line 

11 

H 

1.081797  0 

121.521232 

0 

-60.000000 

0 

1 

2 

3 

Line 

12 

0 

0.000000  0 

0.000000 

0 

0.000000 

0 

0 

0 

0 

Line 

13 

3  1 

4  5  6  7  8 

Line 

14 

3  2 

4  5  6  7  8 

Line 

15 

6  7 

7 

Line 

16 

6  11 

8 

Line 

17 

Line 

18 

70 

80  90  100  no 

120  130  140 

150 

Points  to  note: 

(1)  The  dihedrals  of  the  second  and  third  hydrogens  are  not  marked 
for  optimization:  the  dihedrals  follow  from  point-group  symmetry. 

(2)  All  six  C-H  bond  lengths  and  H-C-C  angles  are  related  by 
symmetry:  see  lines  13  and  14. 

(3)  The  dihedral  on  line  9  Is  the  reaction  coordinate,  while  the 
dihedrals  on  lines  10  and  11  are  related  to  It  by  symmetry  functions  on 
lines  15  and  16  (see  SYMMETRY  for  definitions  of  functions  1,  2,  7,  and 
11). 

(4)  Symmetry  data  are  ended  by  a  blank  line. 

(5)  The  reaction  coordinate  data  are  ended  by  the  end  of  file; 
several  lines  of  data  are  allowed. 


! 
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CHAPTER  6 


BACKGROUND 


6.1  INTRODUCTION 

While  all  the  theory  used  in  MOPAC  Is  In  the  literature,  so  that  In 
principle  one  to  could  read  and  understand  the  algorithm,  many  parts  of 
the  code  Involve  programming  concepts  or  constructions  which,  while  not 
of  sufficient  importance  to  warrant  publication,  are  described  here  In 
order  to  facilitate  understanding. 


6.2  CONVERGENCE  TESTS  IN  SUBROUTINE  ITER 

Self-Consistency  Test 

The  SCF  Iterations  are  stopped  when  two  tests  are  satisfied.  These 
are  (1)  when  the  difference  In  electronic  energy,  In  eV,  between  any  two 
consecutive  Iterations  drops  below  SELCON  and  the  difference  between  any 
three  consecutive  Iterations  drops  below  ten  times  SELCON,  and  (2)  the 
difference  In  density  matrix  elements  on  two  successive  Iterations  falls 
below  a  preset  limit,  which  Is  a  multiple  of  SELCON. 

SELCON  Is  set  initially  to  0.00001  kcal/mole;  this  can  be  made  100 
times  smaller  by  specifying  PRECISE  or  FORCE.  It  can  be  over-ridden  by 
explicitly  defining  the  SCF  criterion  via  SCFCRT=0. 0000001 . 

SELCON  Is  further  modified  by  the  value  of  the  gradient  norm.  If 
known.  If  GN0RM  is  large,  then  a  more  lax  SCF  criterion  is  acceptable, 
and  SCFCRT  can  be  relaxed  up  to  50  times  Its  default  value.  As  the 
gradient  norm  drops,  the  SCF  criterion  returns  to  Its  default  value. 

The  SCF  test  Is  performed  using  the  energy  calculated  from  the  Fock 
matrix  which  arises  from  a  density  matrix,  and  not  from  the  density 
matrix  which  arises  from  a  Fock.  In  the  limit,  the  two  energies  would 
be  Identical,  but  the  first  converges  faster  than  the  second,  without 
loss  of  precision. 
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6.3  CONVERGENCE  IN  SCF  CALCULATION 

A  brief  description  of  the  convergence  techniques  used  In 
subroutine  ITER  follows. 

ITER,  the  SCF  calculation,  employs  six  methods  to  achieve  a 
self-consistent  field.  In  order  of  usage,  these  are: 

(a)  Intrinsic  convergence  by  virtue  of  the  way  the  calculation  Is 
carried  out.  Thus  a  trial  Fock  gives  rise  to  a  trial  density  matrix, 
which  In  turn  Is  used  to  generate  a  better  Fock  matrix. 

This  Is  normally  convergent,  but  many  exceptions  are  known.  The 
main  situations  when  the  Intrinsic  convergence  does  not  work  are: 

(1)  A  bad  starting  density  matrix.  This  normally  occurs 
when  the  default  starting  density  matrix  Is  used.  This  Is  a  very 
crude  approximation,  and  Is  only  used  to  get  the  calculation 
started.  A  large  charge  Is  generated  on  an  atom  In  the  first 
Iteration  the  second  Iteration  overcompensates,  and  an 
oscillation  Is  generated. 

(2)  The  equations  are  only  very  slowly  convergent.  This  can 
be  due  to  a  long-lived  oscillation  or  to  a  slow  transfer  of 
charge. 

(b)  Oscillation  damping.  If,  on  any  two  consecutive  Iterations,  a 
density  matrix  element  changes  by  more  than  0.05,  then  the  density 
matrix  element  Is  set  equal  to  the  old  element  shifted  by  0.05  In  the 
direction  of  the  calculated  element.  Thus,  If  on  Iterations  3  and  4  a 
certain  density  matrix  element  was  0.55  and  0.78,  respectively,  then  the 
element  would  be  set  to  0.60  (-0.55+0.05)  on  Iteration  4.  The  density 
matrix  from  Iteration  4  would  then  be  used  In  the  construction  of  the 
next  Fock  matrix.  The  arrays  which  hold  the  old  density  matrices  are 
not  filled  until  after  Iteration  2.  For  this  reason  they  are  not  used 
In  the  damping  before  Iteration  3. 

(c)  Three-point  Interpolation  of  the  density  matrix.  Subroutine 
CNVG  monitors  the  number  of  Iterations,  and  If  this  Is  exactly  divisible 
by  three,  and  certain  other  conditions  relating  to  the  density  matrices 
are  satisfied,  a  three-point  Interpolation  Is  performed.  This  Is  the 
default  converger,  and  Is  very  effective  with  normally  convergent 
calculations.  It  falls  In  certain  systems,  usually  those  where 
significant  charge  buildup  Is  present. 

(d)  Energy-level  shift  technique.  The  virtual  M.0.  energy  levels 
are  shifted  to  more  positive  energy.  This  has  the  effect  of  damping 
oscillations,  and  Intrinsically  divergent  equations  can  often  be  changed 
to  Intrinsically  convergent  form.  (Invoked  by  the  keyword  SHIFT-n.nnnn) 

(e)  Pulay's  method.  If  requested,  when  the  largest  change  In 
density  matrix  elements  on  two  consecutive  Iterations  has  dropped  below 
0.1,  then  routine  CNVG  Is  abandoned  In  favour  of  a  multl-Fock  matrix 
Interpolation.  This  relies  on  the  fact  that  the  eigenvectors  of  the 
density  and  Fock  matrices  are  Identical  at  self-consistency,  so  [P.F]«0 
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at  SCF.  The  extent  to  which  this  condition  does  not  occur  Is  a  measure 
of  the  deviance  from  self-consistency.  Pul  ay's  method  uses  this 
relationship  to  calculate  that  linear  combination  of  Fock  matrices  which 
minimize  tP.F].  This  new  Fock  matrix  Is  then  used  In  the  SCF 
calculation. 

Under  certain  circumstances,  Pulay's  method  can  cause  very  slow 
convergence,  but  sometimes  It  Is  the  only  way  to  achieve  a 
self-consistent  field.  At  other  times  the  procedure  gives  a  ten-fold 
Increase  In  speed,  so  care  must  be  exercised  In  Its  use.  (Invoked  by 
the  keyword  PULAY) 


6.4  CAUSES  OF  FAILURE  TO  ACHIEVE  A  SCF 

In  a  system  where  a  bl radical  can  form,  such  as  ethane  decomposing 
Into  two  CH3  units,  the  normal  RHF  procedure  can  fail  to  go 
self-consistent.  If  the  system  has  marked  blradlcalold  character,  then 
BIRADICAL  or  UHF  and  TRIPLET  can  often  prove  successful.  These  options 
rely  on  the  assumption  that  two  unpaired  electrons  can  represent  the 
open  shell  part  of  the  wave-function. 

Consider  H-Cl ,  with  the  Interatomic  distance  being  steadily 
Increased.  At  first  the  covalent  bond  will  be  strong,  and  a 
self-consistent  field  Is  readily  obtained.  Gradually  the  bond  will 
become  more  Ionic,  and  eventually  the  charge  on  chlorine  will  become 
very  large.  The  hydrogen,  meanwhile,  will  become  very  electropositive, 
and  there  will  be  an  Increased  energy  advantage  to  any  one  electron  to 
transfer  from  chlorine  to  hydrogen.  If  this  In  fact  occurred,  the 
hydrogen  would  suddenly  become  very  electron-rich  and  would,  on  the  next 
Iteration,  lose  its  extra  electron  to  the  chlorine.  A  sustained 
oscillation  would  then  be  Initiated.  To  prevent  this,  If  BIRADICAL  Is 
specified,  exactly  one  electron  will  end  up  on  hydrogen.  A  similar 
result  can  be  obtained  by  specifying  TRIPLET  in  a  UHF  calculation. 


6.5  TORSION  ANGLE  COHERENCY 

MOPAC  calculations  do  not  distinguish  between  enantiomers, 
consequently  the  sign  of  the  dihedrals  can  be  multiplied  by  -1  and  the 
calculations  will  be  unaffected.  However,  if  chirality  is  Important,  a 
user  should  be  aware  of  the  sign  convention  used. 

The  dihedral  angle  convention  used  in  MOPAC  is  that  defined  by 
Klyne  and  Prelog  in  Experientla  16,  521  (1960).  In  this  convention, 
four  atoms,  AXYB,  with  a  dihedral  angle  of  90  degrees,  will  have  atom  B 
rotated  by  90  degrees  clockwise  relative  to  A  when  X  and  Y  are  lined  up 
In  the  direction  of  sight,  X  being  nearer  to  the  eye.  In  their  words, 
"To  distinguish  between  enantiomeric  types  the  angle  :tau:  Is 
considered  as  positive  when  It  Is  measured  clockwise  from  the  front 
substituent  A  to  the  rear  substituent  B,  and  negative  when  it  Is 
measured  anticlockwise."  The  alternative  convention  was  used  In  all 
earlier  programs.  Including  QCPE  353. 
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6.6  VIBRATIONAL  ANALYSIS 

Analyzing  normal  coordinates  Is  very  tedious.  Users  are  normally 
familiar  with  the  Internal  coordinates  of  the  system  they  are  studying, 
but  not  familiar  with  the  cartesian  coordinates.  To  help  characterize 
the  normal  coordinates,  a  very  simple  analysis  Is  done  automatically, 
and  users  are  strongly  encouraged  to  use  this  analysis  first,  and  then 
to  look  at  the  normal  coordinate  eigenvectors. 

In  the  analysis,  each  pair  of  bonded  atoms  Is  examined  to  see  If 
there  Is  a  large  relative  motion  between  them.  By  bonded  Is  meant 
within  the  Van  der  Haals'  distance.  If  there  Is  such  a  motion,  the 
Indices  of  the  atoms,  the  relative  distance  In  Angstroms,  and  the 
percentage  radial  motion  are  printed.  Radial  plus  tangential  motion 
adds  to  100%,  but  as  there  are  two  orthogonal  tangential  motions  and 
only  one  radial,  the  radial  component  Is  printed. 


6.7  SPARKLES 


Two  extra  "elements"  have  been  put  Into  MOPAC.  These  represent 
pu. e  Ionic  charges,  roughly  equivalent  to  the  following  chemical 
entitles: 

Chemical  Symbol  Equivalent  to 


+ 

atom. 


Tetramethyl  ammonium  radical,  Potassium  atom  or  Cesium 
Borohydrlde  radical.  Halogen,  or  Nitrate  radical 


For  the  purposes  of  discussion  these  entitles  are  called 
"sparkles":  the  name  arises  from  consideration  of  their  behavior. 

Behavior  of  sparkles  In  MOPAC. 

Sparkles  have  the  following  properties: 

1.  Their  nuclear  charge  Is  Integer,  and  Is  +1  or  -1;  there  are  an 
equivalent  number  of  electrons  to  maintain  electroneutral Ity,  1 
and  -1  respectively. 

2.  They  all  have  an  Ionic  radius  of  0.7  Angstroms.  Any  two 
sparkles  of  opposite  sign  will  form  an  Ion-pair  with  a 
Interatomic  separation  of  1.4A. 

3.  They  have  a  zero  heat  of  atomization,  no  orbitals,  and  no 
Ionization  potential . 


They  can  be  regarded  as  unpolarlzable  Ions  of  diameter  1.4A.  They 
do  not  contribute  to  the  orbital  count,  and  cannot  accept  or  donate 
electrons. 
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Since  they  appear  as  uncharged  species  which  Immediately  ionize, 
attention  should  be  given  to  the  charge  on  the  whole  system.  For 
example,  If  the  alkaline  metal  salt  of  formic  acid  was  run,  the  formula 
would  be: 

HC00+  where  +  Is  the  unipositive  sparkle.  The  charge  on  the  system 
would  then  be  zero. 

A  water  molecule  polarized  by  a  positive  sparkle  would  have  the 
formula  H20+,  and  the  charge  on  the  system  would  be  +1 

At  first  sight,  a  sparkle  would  appear  to  be  too  ionic  to  be  a 
point  charge  and  would  combine  with  the  first  charge  of  opposite  sign  it 
encountered. 

This  representation  is  faulty,  and  a  better  description  would  be  of 
an  Ion,  of  diameter  1.4A,  and  the  charge  delocalized  over  Its  surface. 
Computationally,  a  sparkle  Is  an  Integer  charge  at  the  center  of  a 
repulsion  sphere  of  form  exp(-alpha*r) .  The  hardness  of  the  sphere  Is 
such  that  other  atoms  or  sparkles  can  approach  within  about  2  Angstroms 
quite  easily,  but  only  with  great  difficulty  come  closer  than  1.4A. 

Uses  of  Sparkles 

1.  They  can  be  used  as  counterions,  e.g.  for  acid  anions  or  for 
cations.  Thus,  If  the  Ionic  form  of  an  acid  Is  wanted,  then 
the  moieties  H.X,  H.-,  and  +.X  could  be  examined. 

2.  Two  sparkles  of  equal  and  opposite  sign  can  form  a  dipole  for 
mimicking  solvation  effects.  Thus  water  could  be  surrounded  by 
six  dipoles  to  simulate  the  solvent  cage.  A  dipole  of  value  D 
can  be  made  by  using  the  two  sparkles  +  and  -.  If  +  and  -  are 
used,  the  Inter-sparkle  separation  would  be  D/4.803  Angstroms. 
If  the  Inter-sparkle  separation  Is  less  than  1.0  Angstroms  (a 
situation  that  cannot  occur  naturally)  then  the  energy  due  to 
the  dipole  on  Its  own  is  subtracted  from  the  total  energy. 


6.8  MECHANISM  OF  THE  FRAME  IN  THE  FORCE  CALCULATION 

The  FORCE  calculation  uses  cartesian  coordinates,  and  all  3N  modes 
are  calculated,  where  N  Is  the  number  of  atoms  In  the  system.  Clearly, 
there  will  be  5  or  6  "trivial"  vibrations,  which  represent  the  three 
translations  and  two  or  three  rotations.  If  the  molecule  Is  exactly  at 
a  stationary  point,  then  these  "vibrations"  will  have  a  force  constant 
and  frequency  of  precisely  zero.  If  the  force  calculation  was  done 
correctly,  and  the  molecule  was  not  exactly  at  a  stationary  point,  then 
the  three  translations  should  be  exactly  zero,  but  the  rotations  would 
be  non-zero.  The  extent  to  which  the  rotations  are  non-zero  Is  a 
measure  of  the  error  In.  the  geometry. 
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If  the  distortions  are  non-zero,  the  trivial  vibrations  can 
Interact  with  the  low-lying  genuine  vibrations  or  rotations,  and  with 
the  transition  vibration  If  present. 

To  prevent  this  the  analytic  form  of  the  rotations  and  vibrations 
Is  calculated,  and  arbitrary  eigenvalues  assigned;  these  are  500,  600, 
700,  800,  900,  and  1000  ml  1 1 Idynes/angstrom  for  Tx,  Ty,  Tz,  Rx,  Ry  and 
Rz  (If  present),  respectively.  The  rotations  are  about  the  principal 
axes  of  Inertia  for  the  system,  taking  Into  account  Isotopic  masses. 
The  "force  matrix"  for  these  trivial  vibrations  Is  determined,  and  added 
on  to  the  calculated  force  matrix.  After  diagonal Izatlon  the  arbitrary 
eigenvalues  are  subtracted  off  the  trivial  vibrations,  and  the  resulting 
numbers  are  the  "true"  values.  Interference  with  genuine  vibrations  Is 
thus  avoided. 


6.9  PSEUD0DI AG0NALI ZATI0N  -SUBROUTINE  DIAG 

The  basis  of  subroutine  DIAG  Is  the  observation  that  accurate 
matrix  diagonal Izatlon  of  the  secular  determinant  Is  not  a  prerequisite 
In  the  SCF  procedure  for  obtaining  a  self-consistent  density  matrix  In  a 
variational 1y  optimized  calculation.  To  have  a  self-consistent  density 
matrix  It  Is  sufficient  to  have  annihilated  all  energy  matrix  elements 
connecting  the  occupied  and  virtual  molecular  orbitals. 

THEORY 

Given  a  basis  set  of  N  atomic  orbitals  and  Ne  electrons,  there  will 
be  No  -  Ne  /2  occupied  molecular  orbitals  and  Nvlr  «  N-Nv  virtual 
orbitals.  If  the  approximate  form  of  the  molecular  orbitals  Is  known, 
perhaps  from  an  accurate,  standard  diagonal Izatlon  of  the  first  trial 
secular  determinant,  then  the  Interaction  matrix  can  be  constructed. 
The  off-diagonal  matrix  elements  In  Fov  can  then  be  annihilated  by  a 
series  of  2x2  rotations  In  the  manner  of  Jacobi.  Unlike  Jacobi, 
however,  the  method  need  not  be  cycled  to  exactly  diagonalize  Fov;  It  Is 
sufficient  to  have  only  one  sweep.  This  Is  due  to  the  fact  that  the 
Fock  equations  form  a  pseudo-eigenvalue  problem,  and  It  Is  necessary  to 
Iterate  to  obtain  a  self-consistent  field.  For  the  same  reason  the 
second-order  effects  of  the  2x2  rotations  can  be  Ignored.  In  an  exact 
diagonal Izatlon  the  off-diagonal  matrix  elements  formed  by  an  elementary 
2  by  2  rotation  would  have  to  be  eliminated.  These  are  normally  less 
than  one  tenth  of  the  matrix  element  being  annihilated,  and  as  the  SCF 
procedure  does  not  converge  at  one  magnitude  per  cycle  the  second-order 
errors  Introduced  can  be  absorbed  into  the  Fock  matrix  of  the  following 
cycle.  Also,  since  second-order  effects  In  the  "diagonal Izatlon"  are 
being  Ignored  It  Is  equally  valid  to  eliminate  only  those  matrix 
elements  which  are  comparable  with  the  largest  off-diagonal  elements  In 
Fov.  A  further  advantage  of  the  pseudo-eigenvalue  nature  of  the  SCF 
equations  appears  when  we  come  to  evaluate  the  diagonal  terms  of  the 
secular  determinant.  For  this,  we  can  equate  these  elements  with  the 
eigenvalues  resulting  from  the  exact  diagonal Izatlon,  and  hold  them 
exactly  constant  throughout  the  entire  calculation,  right  up  to 
self-consistency.  At  first  sight  this  would  appear  to  introduce  errors 
In  the  final  SCF  density  matrix,  as  obviously  the  sum  of  the  eigenvalues 
cannot  be  constant  In  an  exact  calculation,  and  thus  the  final  sum  of 
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occupied  energy  levels  must  be  in  error.  However,  to  obtain  a  SCF 
density  matrix  we  not  only  do  not  need  to  know  the  exact  eigenvalues,  we 
have  no  need  to  know  the  sum  of  the  occupied  energy  levels.  Using  the 
Initial  set  of  eigenvalues,  a  2  by  2  rotation  will,  of  course,  not 
eliminate  fully  even  those  elements  which  we  do  choose  to  operate  on, 
but  again  the  pseudoeigenvalue  nature  of  the  problem  comes  to  our 
rescue.  As  the  Iterations  proceed  those  errors  introduced  are  rapidly 
eliminated,  so  that  at  self-consistency  an  exact  density  matrix  Is 
generated,  but  we  have  no  knowledge  of  the  values  of  the  eigenvalues, 
eigenvectors  or  two-electron  energy.  This  completes  the  definition  of 
the  secular  determinant,  and  the  2  by  2  rotations  needed  to 
pseudo-diagonal Ize  it.  The  unitary  matrix  that  results  from  the  set  of 
rotations  would  normally  need  to  be  multiplied  by  the  original  set  of 
eigenvectors  to  obtain  the  correct  molecular  orbital  matrix.  There  Is 
no  reason  to  start  the  pseudo-diagonal izatlon  with  a  unit  matrix,  and 
this  step  can  be  eliminated  by  starting  the  pseudo-diagonal Izatlon  using 
the  old  set  of  molecular  orbitals. 
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COMPUTATIONAL  ASPECTS 

Clearly,  an  Initial  exact  diagonal Izatlon  Is  needed  In  order  to 
obtain  a  good  starting  set  of  eigenvectors.  This  Is,  however,  not  a 
sufficient  condition  for  Initiating  the  use  of  the  new  method.  During 
the  first  one  or  two  cycles  of  a  SCF  procedure  the  order  of  occupancy  of 
the  molecular  orbitals  may  change.  That  Is,  the  occupancy  of  the  M.O.'s 
In  the  first  two  Iterations  may  correspond  to  an  excited  singlet  state. 
For  this  reason  It  Is  recommended  that  the  new  method  be  used  only  after 
the  Initial  large  fluctuations  In  the  density  matrix  have  died  down. 
Our  arbitrary  criterion  was  that  the  largest  change  In  the  diagonal 
elements  of  the  density  matrix  between  any  two  Iterations  should  be  less 
than  0.05. 

The  method  proposed  uses  three  matrices:  the  C  matrix  of 
eigenvector  coefficients,  the  F  matrix,  and  a  working  matrix  to  hold  the 
new  Fov  matrix.  In  virtual  memory  computers  this  would  normally  not 
present  any  problem.  A  criterion  Is  also  needed  to  decide  when  to 
perform  a  2  by  2  rotation  and  elimination.  After  a  few  trial 
calculations.  It  was  found  to  be  efficient  to  eliminate  all  off-diagonal 
elements  whose  modulus  was  larger  than  0.01  times  the  modulus  of  the 
largest  off-diagonal  element.  No  Important  change  In  the  number  of 
Iterations  was  observed.  The  electronic  energy  calculated  after  each 
Iteration  does  Indeed  differ  following  the  Introduction  of  the  new 
method.  This  Is  a  result  of  all  the  second-order  effects  being 
Introduced,  but  the  differences  between  electronic  energies  calculated 
by  the  exact  diagonal Izatlon  and  by  the  new  method  rapidly  converge.  An 
Interesting  side-effect  of  the  new  non-rigorous  method  Is  that  some 
damping  Is  Introduced:  In  all  molecules  examined  the  number  of 
Iterations  required  to  achieve  self-consistency  either  stayed  constant 
or  dropped  by  one  or  two. 
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6.10  CONFIGURATION  INTERACTION 

MOPAC  contains  a  set  of  subroutines  capable  of  a  very  large 
Multi-Electron  Configuration  Interaction  calculation,  MECI,  which  allows 
almost  any  configuration  Interaction  calculation  to  be  performed. 
Because  of  Its  complexity,  two  distinct  levels  of  input  are  supported; 
the  default  values  will  be  of  use  to  the  novice  while  an  expert  has 
available  an  exhaustive  set  of  key-words  from  which  a  specific  C.I.  can 
be  tal lored. 

A  MECI  calculation  involves  the  interaction  of  microstates 
representing  specific  permutations  of  electrons  in  a  set  of  M.O.s. 
Starting  with  a  set  electronic  configuration,  either  closed  shell  or 
open  shell,  but  unconditionally  restricted  Hartree-Fock,  the  first  step 
In  a  MECI  calculation  Is  the  removal  from  the  M.O.s  of  the  electrons  to 
be  used  In  the  C.I. 

Each  microstate  is  then  constructed  from  these  empty  M.O.s  by 
adding  In  electrons  according  to  a  prescription.  The  energy  of  the 
configuration  Is  evaluated,  as  is  the  energy  of  interaction  with  all 
previously-defined  configurations.  Diagonal Izatlon  then  results  In 
state  functions.  From  the  eigenvectors  the  expectation  value  of  s**2  Is 
calculated,  and  the  spin-states  of  the  state  functions  calculated. 


General  Overview  of  Key-Words 


Key-words  associated  with  the  operations  of  MECI  are: 
SINGLET  DOUBLET  EXCITED 

TRIPLET  BIRADICAL  ESR 

OPEN(nl,n2)  C.I. =2  MECI 

R00T=n 


Each  key-word  may  Imply  others;  thus  TRIPLET  Implies  an  open-shell 
system,  therefore  OPEN<2,2),  and  C.I. =2  are  implied.  If  not  user 
specified. 


Starting  Electronic  Configuration 


MECI  Is  restricted  to  RHF  calculations,  but  with  that  single 
restriction  any  starting  configuration  will  be  supported.  Examples  of 
starting  configurations  would  be 

System  Key-Words  used  Starting  Configuration 


Methane 
Methyl  Radical 
Twisted  Ethylene 
Twisted  Ethylene 
Twisted  Ethylene 


<none> 
<none> 
TRIPLET 
•0PEN(2 ,2 ) 
Cation  OPENd  ,2) 


2.00  2.00  2.00  2.00  2.00 
2.00  2.00  2.00  2.00  1 .00 
2.00  2.00  2.00  1 .00  1 .00 
2.00  2.00  2.00  1 .00  1.00 
2.00  2.00  2.00  0.50  0.50 
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Choice  of  starting  configuration  Is  Important.  For  example.  If 
twisted  ethylene,  a  ground-state  triplet,  Is  not  defined  using  TRIPLET 
or  0PEN(2,2),  then  the  closed-shell  ground-state  structure  will  be 
calculated.  Obviously  this  configuration  Is  a  legitimate  microstate, 
but  from  the  symmetry  of  the  system  a  better  choice  would  be  to  define 
one  electron  In  each  of  the  two  formally  degenerate  pi-type  M.O.s.  The 
Initial  SCF  calculation  does  not  distinguish  between  OPEN(2,2)  and 
TRIPLET  since  both  key-words  define  the  same  starting  configuration. 
This  can  be  verified  by  monitoring  the  convergance  using  PL,  for  which 
both  key-words  give  the  same  SCF  energy. 


Removal  of  Electrons  from  Starting  Configuration 

For  a  starting  configuration  of  alpha  M.O.  occupancies  0(1),  0(1) 

being  In  the  range  0.0  to  1.0,  the  energies  of  the  M.O.s  Involved  In  the 
MECI  can  be  calculated  from 

E(1)  -  Sum(j)(2J(1 ,j)-K(1 ,j))0(j) 


where  JO, 3)  and  K(1,j)  are  the  coulomb  and  exchange  Integrals  between 
M.O.s  1  and  j.  The  M.O.  Index  j  runs  over  those  M.O.s  Involved  In  the 
MECI  only.  Most  MECI  calculations  will  Involve  between  1  and  5  M.O.s, 
so  a  system  with  about  30  filled  or  partly  filled  M.O.s  could  have  M.O.s 
25-30  Involved.  The  resulting  eigenvalues  correspond  to  those  of  the 
cationic  system  resulting  from  removal  of  n  electrons,  where  n  Is  twice 
the  sum  of  the  orbital  occupancies  of  those  M.O.s  Involved  In  the  C.I. 

The  arbltary  zero  of  energy  In  a  MECI  calculation  Is  the  starting 
ground  state,  without  any  correction  for  errors  Introduced  by  the  use  of 
fractional  occupancies.  In  order  to  calculate  the  energy  of  the  various 
configurations,  the  energy  of  the  vacuum  state  ( 1 . e . ,  the  state 
resulting  from  removal  of  the  electrons  used  In  the  C.I.)  needs  to  be 
evaluated.  This  energy  Is  defined  by 


GSE  -  Sum(1)[  E(  1  )0(  1 )  +  J(1,1>  *  0(1  >*0(1  > 

+  Sum( j<1 )(  2(2J( 1  )  -  K(1 ,j )  )  *  0(1 )*0(3>  >  ] 


Multi  Electron  Configuration  Interaction 


The  numbering  of  the  M.O.s  used  In  the  MECI  is  standard,  and 
follows  the  Aufbau  principle.  The  order  of  filling  Is  In  order  of 
energy,  and  alpha  before  beta.  This  pofnt  Is  critically  important  in 
deciding  the  sign  of  matrix  elements.  For  a  5  M.O.  system,  then,  the 
order  of  filling  Is. 
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(1X1X2X2X3X3X4X4X5X5) 


A  triplet  state  arising  from  two  microstates,  each  with  a  component 
of  spin  -  0,  will  thus  be  the  positive  combination. 

<I)<2)  +  (1X2) 


This  Is  In  variance  with  the  sign  convention  used  In  earlier 
programs  for  running  MNDO.  This  standard  sign  convention  was  chosen  In 
order  to  allow  the  signs  of  the  microstate  coefficients  to  conform  to 
those  resulting  from  the  spin  step-down  operator. 

Matrix  elements  between  all  pairs  of  microstates  are  calculated  in 
order  to  form  the  secular  determinant.  Many  elements  will  be 
Identically  zero,  due  to  the  Interacting  determinants  differing  by  more 
than  two  M.O.s.  For  the  remaining  Interactions  the  following  types  can 
be  Identified. 

1.  The  two  determinants  are  identical: 

No  permutations  are  necessary  In  order  to  calculate  the 
sign  of  the  matrix  element.  E(p,p)  Is  given  simply  by 

E(p,p)«Sum(1)0a(1 ,p)*[E1g(1)  +  Sum( 1/2(<1 1£j j>-<1 j£1j>)*0a< j ,p> 

+■  <11ijJ>  *0b(j,p>) 

+  Sum( 1 )0b( 1 ,p)*[E1g( 1 )  ♦  Sum(l/2(<1 Iijj>-<1 j£1J>)*Ob( j ,p>> 


0a(1,p)  «  Occupancy  of  alpha  M.O.  i  In  Microstate  p 
0b(1,p)  -  Occupancy  of  beta  M.O.  1  In  Microstate  p 


2.  Determinants  differing  by  exactly  one  M.O.: 

The  differing  M.O.  can  be  of  type  alpha  or  beta.  It  Is 
.  sufficient  to  evaluate  the  case  In  which  both  M.O.s  are  of 
alpha  type,  the  beta  form  Is  obtained  In  like  manner. 

E(p,q)  -  Sum(k)  [  <1j£kk>  -  <1k£jk>  )  *  (Occa(k)  -  Occg(k)) 

+  <1jikk>  *  (Occb(k)  -  Occg(k) ] 

E(p,q)  may  need  to  be  multiplied  by  -1,  If  the  number  of 

two  electron  permutations  required  to  bring  M.O.s  1  and  j 
Into  coincidence  Is  odd. 


Where  Occa(k)  Is  the  alpha  molecular  orbital  occupancy  In 
the  configuration  Interaction. 
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3.  Determinants  differing  by  exactly  two  M.O.s: 

The  two  M.O.s  can  have  the  same  or  opposite  spins.  Three 

cases  can  be  Identified: 

1.  Both  M.O.s  have  alpha  spin: 

For  the  first  microstate  having  M.O.s  1  and  J ,  and  the 
second  microstate  having  M.O.s  k  and  1,  the  matrix  element 
connecting  the  two  microstates  is  given  by 

Q(p,q)  *  <1k£jl >  -  <1 l£jk> 

E(p,q)  may  need  to  be  multiplied  by  -1,  If  the  number  of 
two  electron  permutations  required  to  bring  M.O.  1  Into 
coincidence  with  M.O.  k  and  M.O.  j  Into  coincidence  with 
M.O.  1  Is  odd. 

2.  Both  M.O.s  have  beta  spin: 

The  matrix  element  is  calculated  In  the  same  manner  as 
In  the  previous  case. 

3.  One  M.O.  has  alpha  spin,  and  one  beta  spin: 

For  the  first  microstate  having  M.O.s  alpha(l)  and 
beta(j),  and  the  second  microstate  having  M.O.s  alpha(k) 
and  beta(l),  the  matrix  element  connecting  the  two 
microstates  Is  given  by 

Q(p,q)  =  <1k£jl > 

E(p,q)  may  need  to  be  multiplied  by  -1,  If  the  number  of 
two  electron  permutations  required  to  bring  M.O.  i  into 
coincidence  with  M.O.  k  and  M.O.  j  into  coincidence  with 
M.O.  1  Is  odd. 


States  Arising  from  Various  Calculations 


Each  MECI  calculation  Invoked  by  use  of  the  key-word  C.I.-n 
normally  gives  rise  to  states  of  quantized  spins.  When  C.I.  Is  used 
without  any  other  modifying  key-words,  the  following  states  will  be 
obtained. 
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No.  of  M.O.s  States  Arising  States  Arising  From 

From  Odd  Electron  Systems  Even  Electron  Systems 
in  MECI  Doublets  Singlets  Triplets 

1  1  1 

2  2  3  1 

These  numbers  of  spin  states  will  be  obtained  irrespective  of  the 
chemical  nature  of  the  system. 


Calculation  of  Spin-States 


In  order  to  calculate  the  spin-state,  the  expectation  value  of  S2 
is  calculated. 

S2  =  SCS+1)  =  Sz**2  +  2*S(+)S(-) 

=  Ne  - 

Sum(1>  [CO  ,k)*C( i  ,k)*(l  /4*(Na(  i  )-Nb(  1 )  >**2 
+  Sum(l)  Oa(l ,i )*Ob(l ,i )) 

+Sum<j)  2[C<1,k>*C<j,k)*<Kronekerdelta(C(i,k>(  S(+)S(->  )C(j,k)]] 

Where  Ne  =  No.  of  electrons  in  C.I. 

CO.k)  *  Coefficient  of  Microstate  1  in  State  k 

Na(i)  =  Number  of  alpha  electrons  in  Microstate  i 

Nb(1)  *  Number  of  beta  electrons  in  Microstate  1 

Oa(l,k)  *  Occupancy  of  alpha  M.O.  1  in  Microstate  k 

Ob(l.k)  =  Occupancy  of  beta  M.O.  1  in  Microstate  k 

S(+)  =  Spin  shift  up  or  step  up  operator 

S(-)  =  Spin  shift  down  or  step  down  operator 

The  Kronekerdel ta  Is  1  if  the  two  terms  in  brackets  following  it 

are  Identical. 


The  spin  state  is  calculated  from  S  =  1/2  (  Sqrt(l+4*S2>  -  1  ) 

In  practice,  S  is  calculated  to  be  exactly  Integer,  or  half 
Integer.  That  is,  there  is  insignificant  error  due  to  approximations 
used.  This  does  not  mean,  however,  that  the  method  is  accurate.  The 
spin  calculation  is  completely  precise,  in  the  group  theoretic  sense, 
but  the  accuracy  of  the  calculation  is  limited  by  the  Hamiltonian  used, 
a  space-dependent  function. 


Choice  of  State  to  be  Optimized 

MECI  can  calculate  four  states  of  various  total  spin.  Two  schemes 
are  provided  to  allow  a- given  state  to  be  selected.  First,  R00T=n  will, 
when  used  on  its  own,  select  the  n'th  state.  Irrespective  of  its  total 
spin.  By  default  n-1.  If  ROOT-n  Is  used  In  conjunction  with  a  key-word 
from  the  set  SINGLET,  DOUBLET,  or  TRIPLET,  then  the  n'th  root  of  that 
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spin-state  will  be  used.  For  example,  R00T«2  and  SINGLET  will  select 
the  second  singlet  state.  If  there  are  two  triplet  states  below  the 
second  singlet  state  then  this  will  mean  that  the  fourth  state  will  be 
selected. 


Calculation  of  Unpaired  Spin  Density 

Starting  with  the  state  functions  as  linear  combinations  of 
configurations,  the  unpaired  spin  density,  corresponding  to  the  alpha 
spin  density  minus  the  beta  spin  density,  will  be  calculated  for  the 
first  few  states.  This  calculation  Is  straightforward  for  diagonal 
terms,  and  only  those  terms  are  used. 
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6.11  REDUCED  MASSES  IN  A  FORCE  CALCULATION 

Reduced  masses  for  a  diatomic  are  given  by 
(massl)  *  (mass2> 


(mass!)  +  <mass2> 

For  a  Hydrogen  molecule  the  reduced  mass  Is  thus  0.5;  for  heavily 
hydrogenated  systems,  e.g.  methane,  the  reduced  mass  can  be  very  low. 
A  vibration  Involving  only  heavy  atoms  ,  e.g.  a  C-N  In  cyanide,  should 
give  a  large  reduced  mass. 

For  the  "trivial"  vibrations  the  reduced  mass  Is  Ill-defined,  and 
where  this  happens  the  reduced  mass  Is  set  to  zero. 


6.12  USE  OF  SADDLE  CALCULATION 

A  SADDLE  calculation  uses  two  complete  geometries,  as  shown  on  the 
following  data  file  for  the  ethyl  radical  hydrogen  migration  from  one 
methyl  group  to  the  other. 
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Details  of  the  mathematics  of  SADDLE  appeared  In  print  In  1984,  so 
only  a  superficial  description  will  be  given  here. 

The  main  steps  In  the  saddle  calculation  are  as  follows: 

1.  The  heats  of  formation  of  both  systems  are  calculated. 

2.  A  vector  R  of  length  3N-6  defining  the  difference  between  the 
two  geometries  Is  calculated. 
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3.  The  scalar  P  of  the  difference  vector  Is  reduced  by  some 
fraction,  normally  about  5  to  15  percent. 

4.  Identify  the  geometry  of  lower  energy;  call  this  G. 

5.  Optimize  G,  subject  to  the  constraint  that  It  maintains  a 
constant  distance  P  from  the  other  geometry. 

6.  If  the  newly-optimized  geometry  Is  higher  In  energy  then  the 
other  geometry,  then  go  to  1 .  If  It  Is  higher,  and  the  last 
two  steps  Involved  the  same  geometry  moving,  make  the  other 
geometry  G  without  modifying  P,  and  go  to  5. 

7.  Otherwise  go  back  to  2. 


The  mechanism  of  5  Involves  the  coordinates  of  the  moving  geometry 
being  perturbed  by  an  amount  equal  to  the  product  of  the  discrepancy 
between  the  calculated  and  required  P  and  the  vector  R. 

As  the  specification  of  the  geometries  Is  quite  difficult,  in  that 
the  difference  vector  depends  on  angles  (which  are,  of  necessity 
Ill-defined  by  360  degrees)  SADDLE  can  be  made  to  run  In  cartesian 
coordinates  using  the  key-word  XYZ.  If  this  option  Is  chosen  then  the 
Initial  steps  of  the  calculation  are  as  follows: 

1.  Both  geometries  are  converted  Into  cartesian  coordinates. 

2.  Both  geometries  are  centered  about  the  origin  of  cartesian 

space. 

3.  One  geometry  Is  rotated  until  the  difference  vector  Is  a 

minimum  -  this  minimum  Is  within  1  degree  of  the  absolute 
bottom. 

4.  The  SADDLE  calculation  then  proceeds  as  described  above. 
LIMITATIONS:  . 

The  two  geometries  must  be  related  by  a  continuous  deformation  of 
the  coordinates.  By  default.  Internal  coordinates  are  used  In 

specifying  geometries,  and  while  bond  lengths  and  bond  angles  are 
unambiguously  defined  (being  both  positive),  the  dihedral  angles  can  be 
positive  or  negative.  Clearly  300  degrees  could  equally  well  be 
specified  as  -60  degrees.  A  wrong  choice  of  dihedral  would  mean  that 
Instead  of  the  desired  reaction  vector  being  used,  a  completely 
Incorrect  vector  was  used,  with  disastrous  results. 

To  correct  this,  ensure  that  one  geometry  can  be  obtained  from  the 
other  by  a  continuous  deformation,  or  use  the  XYZ  option. 


CHAPTER  7 


PROGRAM 


The  logic  within  MOPAC  Is  best  understood  by  use  of  flow-diagrams. 

There  are  two  main  sequences,  geometric  and  electronic.  These  join 
only  at  one  common  subroutine  COMPFG.  It  Is  possible,  therefore,  to 
understand  the  geometric  or  electronic  sections  In  Isolation,  without 
having  studied  the  other  section. 
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7.1  MAIN  GEOMETRIC  SEQUENCE 


COMPFG  (See  ELECTRONIC  SEQUENCE) 


PROGRAM 


7.2  MAIN  ELECTRONIC  FLOW 
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7.3  CONTROL  WITHIN  MOPAC 

Almost  all  the  control  Information  Is  passed  via  the  single  datum 
"KEYWRD",  a  string  of  80  characters,  which  Is  read  In  at  the  start  of 
the  Job. 

Each  subroutine  Is  made  Independent,  as  far  as  possible,  even  at 
the  expense  of  extra  code  or  calculation.  Thus,  for  example,  the  SCF 
criterion  Is  set  In  subroutine  ITER,  and  nowhere  else.  Similarly 
subroutine  DERIV  has  exclusive  control  of  the  step  size  In  the 
finite-difference  calculation  of  the  energy  derivatives.  If  the  default 
values  are  to  be  reset,  then  the  new  value  Is  supplied  In  KEYWRD,  and 
extracted  via  INDEX  and  READA.  The  flow  of  control  Is  decided  by  the 
presence  of  various  key-words  In  KEYWRD. 

When  a  subroutine  Is  called.  It  assumes  that  all  data  required  for 
Its  operation  are  available  In  either  common  blocks  or  arguments. 
Normally  no  check  Is  made  as  to  the  validity  of  the  data  received.  All 
data  are  "owned"  by  one,  and  only  one,  subroutine.  By  ownership  Is 
Implied  the  permission  and  ability  to  change  the  data.  Thus  MOLDAT 
"owns"  the  number  of  atomic  orbitals.  In  that  It  calculates  this  number, 
and  stores  It  In  the  variable  NORBS.  Many  subroutines  use  NORBS,  but 
none  of  them  Is  allowed  to  change  It.  For  obvious  reasons  no  exceptions 
should  be  made  to  this  rule.  To  Illustrate  the  usefulness  of  this 
convention,  consider  the  eigenvectors,  C  and  CBETA.  These  are  owned  by 
ITER.  Before  ITER  Is  called,  C  and  CBETA  are  not  calculated,  after  ITER 
has  been  called  C  and  CBETA  are  known,  so  any  subroutine  which  needs  to 
use  the  eigenvectors  can  do  so  In  the  certain  knowledge  that  they  exist. 

Any  variables  which  are  only  used  within  a  subroutine  are  not 
passed  outside  the  subroutine  unless  an  overriding  reason  exists.  This 
Is  found  In  PULAY  and  CNVG,  among  others  where  arrays  used  to  hold 
spin-dependent  data  are  used,  and  these  cannot  conveniently  be  defined 
within  the  subroutines.  In  these  examples,  the  relevant  arrays  are 
"owned"  by  ITER. 

A  general  subroutine,  of  which  ITER  Is  a  good  example,  handles 
three  kinds  of  data:  First,  data  which  the  subroutine  Is  going  to  work 
on,  for  example  the  one  and  two  electron  matrices;  second,  data 
necessary  to  manipulate  the  first  set  of  data,  such  as  the  number  of 
atomic  orbitals;  third,  the  calculated  quantities,  here  the  electronic 
energy,  and  the  density  and  Fock  matrices. 

Reference  data  are  entered  Into  a  subroutine  by  way  of  the  common 
blocks.  This  Is  to  emphasize  their  peripheral  role.  Thus  the  number  of 
orbitals,  while  essential  to  ITER,  Is  not  central  to  the  task  It  has  to 
perform,  and  Is  passed  through  a  common  block. 

Data  the  subroutine  Is  going  to  work  on  are  passed  via  the  argument 
list.  Thus  the  one  and  two  electron  matrices,  which  are  the  main  reason 
for  ITER's  existence,  are  entered  as  two  of  the  four  arguments.  As  ITER 
does  not  own  these  matrices  It  can  use  them  but  may  not  change  their 
contents.  The  other  arguments  are  EE,  the  electronic  energy,  and 
FULSCF,  a  logical.  EE  Is  owned  by  ITER  even  though  It  first  appears 
before  ITER  Is  called.  FULSCF,  on  the  other  hand.  Is  not  owned  by  ITER, 
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and  Is  used,  but  not  changed. 

Sometimes  common  block  data  should  more  correctly  appear  In  an 
argument  list.  This  Is  usually  not  done  In  order  to  prevent  obscuring 
the  main  role  the  subroutine  has  to  perform.  Thus  ITER  calculates  the 
density  and  Fock  matrices,  but  these  are  not  represented  In  the  argument 
list  as  the  calling  subroutine  never  needs  to  know  them;  Instead,  they 
are  stored  In  common. 

SUBROUTINE  GMETRY:  Description  for  programmers. 

GMETRY  has  two  arguments,  GEO  and  COORD.  On  Input  GEO  contains 
either  (a)  Internal  coordinates  or  (b)  cartesian  coordinates.  On  exit 
COORD  contains  the  cartesian  coordinates. 

The  normal  mode  of  usage  Is  to  supply  the  Internal  coordinates.  In 
which  case  the  connectivity  relations  are  found  In  common  block  GEOKST. 

If  the  contents  of  NA( 1 )  Is  zero,  as  required  for  any  normal 
system,  then  the  normal  Internal  to  cartesian  conversion  Is  carried  out. 

If  the  contents  of  NA(1)  Is  99,  then  the  coordinates  found  In  GEO 
are  assumed  to  be  cartesian,  and  no  conversion  Is  made.  This  Is  the 
situation  In  a  FORCE  calculation. 

A  further  option  exists  within  the  Internal  to  cartesian 
conversion.  If  STEP,  stored  In  common  block  REACTN,  Is  non-zero,  then  a 
reaction  path  Is  assumed,  and  the  Internal  coordinates  are  adjusted 
radially  In  order  that  the  "distance"  In  Internal  coordinate  space  from 
the  geometry  specified  In  GEO  Is  STEPP  away  from  the  geometry  stored  In 
GEOA,  stored  In  REACTN. 

During  the  Internal  to  cartesian  conversion,  the  angle  between  the 
three  atoms  used  In  defining  a  fourth  atom  Is  checked  to  ensure  that  It 
Is  not  near  to  0  or  180  degrees.  If  It  Is  near  to  these  angles,  then 
there  Is  a  high  probability  that  a  faulty  geometry  will  be  generated  and 
to  prevent  this  the  calculation  Is  stopped  and  an  error  message  printed. 

NOTE  1:  If  the  angle  Is  exactly  0  or  180  degrees,  then  the 
calculation  Is  not  terminated:  This  Is  the  normal  situation  In  a 
high-symmetry  molecule  such  as  propyne. 

NOTE  2:  The  check  Is  only  made  If  the  fourth  atom  has  a  bond  angle 
which  Is  not  zero  or  180  degrees. 


CHAPTER  8 


ERROR  MESSAGES  PRODUCED  BY  MOPAC 


MOPAC  produces  several  hundred  messages,  all  of  which  are  Intended 
to  be  self-explanatory.  However,  when  an  error  occurs  It  Is  useful  to 
have  more  Information  than  is  given  In  the  standard  messages. 

The  following  alphabetical  list  gives  more  complete  definitions  of 
the  messages  printed. 

AN  UNOPTIMIZABLE  GEOMETRIC  PARAMETER.... 

When  Internal  coordinates  are  supplied,  six  coordinates  cannot  be 
optimized.  These  are  the  three  coordinates  of  atom  1,  the  angle  and 
dihedral  on  atom  2  and  the  dihedral  on  atom  3.  An  attempt  has  been  made 
to  optimize  one  of  these.  This  Is  usually  Indicative  of  a  typographic 
error,  but  might  simply  be  an  oversight.  Either  way,  the  error  will  be 
corrected  and  the  calculation  will  not  be  stopped  here. 

ATOM  NUMBER  nn  IS  ILLDEFINED 

The  rules  for  definition  of  atom  connectivity  are: 

1.  Atom  2  must  be  connected  to  atom  1  (default  -  no  override) 

2.  Atom  3  must  be  connected  to  atom  1  or  2,  and  make  an  angle  with 
2  or  1 . 

3.  All  other  atoms  must  be  defined  In  terms  of  already-defined 

atoms:  these  atoms  must  all  be  different.  Thus  atom  9  might 

be  connected  to  atom  5,  make  an  angle  with  atom  6,  and  have  a 
dihedral  with  atom  7.  If  the  dihedral  was  with  atom  5,  then 
the  geometry  definition  would  be  faulty. 


If  any  of  these  rules  Is  broken,  a  fatal  error  message  Is  printed, 
and  the  calculation  stopped. 
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ATOMIC  NUMBER  nn  IS  NOT  AVAILABLE  ... 

An  element  has  been  used  for  which  parameters  are  not  available. 
Only  If  a  typographic  error  has  been  made  can  this  be  rectified.  This 
check  Is  not  exhaustive,  in  that  even  If  the  elements  are  acceptable 
there  are  some  combinations  of  elements  within  MINDO/3  that  are  not 
allowed.  This  Is  a  fatal  error  message. 

ATOMIC  NUMBER  OF  nn  ? 

An  atom  has  been  specified  with  a  negative  or  zero  atomic  number. 
This  Is  normally  caused  by  forgetting  to  specify  an  atomic  number  or 
symbol.  This  is  a  fatal  error  message. 

ATOMS  nn  AND  nn  ARE  SEPARATED  BY  nn.nnnn  ANGSTROMS. 

Two  genuine  atoms  (not  dummies)  are  separated  by  a  very  small 
distance.  This  can  occur  when  a  complicated  geometry  Is  being 
optimized,  In  which  case  the  user  may  wish  to  continue.  This  can  be 
done  by  using  the  key-word  GEO-OK.  More  often,  however,  this  message 
Indicates  a  mistake,  and  the  calculation  Is,  by  default,  stopped. 

ATTEMPT  TO  GO  DOWNHILL  IS  UNSUCCESSFUL... 

A  quite  rare  message,  produced  by  Bartel's  gradient  norm 
minimization.  Bartel's  method  attempts  to  minimize  the  gradient  norm  by 
searching  the  gradient  space  for  a  minimum.  Apparently  a  minimum  has 
been  found,  but  not  recognized  as  such.  The  program  has  searched  In  all 
(3N-6)  directions,  and  found  no  way  down,  but  the  criteria  for  a  minimum 
have  not  been  satisfied.  No  advice  Is  available  for  getting  round  this 
error. 


BOTH  SYSTEMS  ARE  ON  THE  SAME  SIDE.. 

A  non-fatal  message,  but  still  cause  for  concern.  During  a  SADDLE 
calculation  the  two  geometries  Involved  are  on  opposite  sides  of  the 
transition  state.  This  situation  Is  verified  at  every  point  by 
calculating  the  cosine  of  the  angle  between  the  two  gradient  vectors. 
For  as  long  as  it  Is  negative,  then  the  two  geometries  are  on  opposite 
sides  of  the  T/S.  If,  however,  the  cosine  becomes  positive,  then  the 
assumption  Is  made  that  one  moiety  has  fallen  over  the  T/S  and  Is  now 
below  the  other  geometry.  That  Is,  It  Is  now  further  from  the  T/S  than 
the  other,  temporarily  fixed,  geometry.  To  correct  this.  Identify 
geometries  corresponding  to  points  on  each  side  of  the  T/S.  (Two 
geometries  on  the  output  separated  by  the  message  "SWAPPING...")  and 
make  up  a  new  data-file  using  these  geometries.  This  corresponds  to 
points  on  the  reaction  path  near  to  the  T/S.  Run  a  new  job  using  these 
two  geometries,  but  with  BAR  set  to  a  third  or  a  quarter  of  its  original 
value,  e.g.  8AR-0.05.  This  normally  allows  the  T/S  to  be  located. 
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C.I.  NOT  ALLOWED  WITH  UHF 

There  Is  no  UHF  configuration  Interaction  calculation  In  MOPAC. 
Either  remove  the  key-word  that  Implies  C.I.  or  the  word  UHF. 

CALCULATION  ABANDONED  AT  THIS  POINT 

A  particularly  annoying  message!  In  order  to  define  an  atom's 
position,  the  three  atoms  used  In  the  connectivity  table  must  not 
accidentally  fall  into  a  straight  line.  This  can  happen  during  a 
geometry  optimization  or  gradient  minimization.  If  they  do,  and  If  the 
angle  made  by  the  atom  being  defined  Is  not  zero  or  180  degrees,  then 
Its  position  becomes  Ill-defined.  This  Is  not  desirable,  and  the 
calculation  will  stop  In  order  to  allow  corrective  action  to  be  taken. 
Note  that  If  the  three  atoms  are  In  an  exactly  straight  line,  this 
message  will  not  be  triggered.  The  good  news  Is  that  the  criterion  used 
to  trigger  this  message  was  set  too  coarsely.  The  criterion  has  been 
tightened  so  that  this  message  now  does  not  often  appear.  Geometric 
Integrity  does  not  appear  to  be  compromized. 

CARTESIAN  COORDINATES  READ  IN,  AND  CALCULATION... 

If  cartesian  coordinates  are  read  In,  but  the  calculation  Is  to  be 
carried  out  using  Internal  coordinates,  then  either  all  possible 
geometric  variables  must  be  optimized,  or  none  can  be  optimized.  If 
only  some  are  marked  for  optimization  then  ambiguity  exists.  For 
example.  If  the  "X"  coordinate  of  atom  6  Is  marked  for  optimization,  but 
the  "Y"  Is  not,  then  when  the  conversion  to  Internal  coordinates  takes 
place,  the  first  coordinate  becomes  a  bond-length,  and  the  second  an 
angle.  These  bear  no  relationship  to  the  "X"  or  “Y"  coordinates.  This 
Is  a  fatal  error. 

CARTESIAN  COORDINATES  READ  IN,  AND  SYMMETRY... 

If  cartesian  coordinates  are  read  In,  but  the  calculation  Is  to  be 
carried  out  using  Internal  coordinates,  then  any  symmetry  relationships 
between  the  cartesian  coordinates  will  not  be  reflected  In  the  Internal 
coordinates.  For  example,  If  the  "Y"  coordinates  of  atoms  5  and  6  are 
equal.  It  does  not  follow  that  the  Internal  coordinate  angles  these 
atoms  make  are  equal.  This  Is  a  fatal  error. 

ERROR  DURING  READ  AT  ATOM  NUMBER  .... 

Something  Is  wrong  with  the  geometry  data.  In  order  to  help  find 
the  error,  the  geometry  already  read  In  Is  printed.  The  error  lies 
either  on  the  last  line  of  the  geometry  printed,  or  on  the  next 
(unprinted)  line.  This  Is  a  fatal  error. 
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<<<< - ****  FAILED  TO  ACHIEVE  SCF.  **** - >>>> 

The  SCF  calculation  failed  to  go  to  completion;  an  unwanted  and 
depressing  message  that  unfortunately  appears  every  so  often. 

To  date  three  unconditional  convergers  have  appeared  In  the 
literature:  the  SHIFT  technique,  Pulay's  method,  and  the  Camp-King 
converger.  It  would  not  be  fair  to  the  authors  to  condemn  their 
methods.  In  MOPAC  all  sorts  of  weird  and  wonderful  systems  are 
calculated,  systems  the  authors  of  the  convergers  never  dreamed  of. 
MOPAC  uses  a  combination  of  all  three  convergers  at  times.  Normally 
only  a  quadratic  damper  Is  used. 

If  this  message  appears,  suspect  first  that  the  calculation  might 
be  faulty,  then.  If  you  feel  confident,  try  altering  the  SHIFT,  or 
Involklng  PULAY  or  CAMP-KING  on  their  own. 

If  nothing  works,  then  consider  slackening  the  SCF  criterion.  This 
will  allow  heats  of  formation  to  be  calculated  with  reasonable 
precision,  but  the  gradients  are  likely  to  be  Imprecise. 

GEOMETRY  TOO  UNSTABLE  FOR  EXTRAPOLATION.. 

In  a  reaction  path  calculation  the  Initial  geometry  for  a  point  Is 
calculated  by  quadratic  extrapolation  using  the  previous  three  points. 

If  a  quadratic  fit  Is  likely  to  lead  to  an  Inferior  geometry,  then 
the  geometry  of  the  last  point  calculated  will  be  used.  The  total 
effect  Is  to  slow  down  the  calculation,  but  no  user  action  Is 
recommended. 

**  GRADIENT  IS  TOO  LARGE  TO  ALLOW... 

Before  a  FORCE  calculation  can  be  performed  the  gradient  norm  must 
be  so  small  that  the  third  and  higher  order  components  of  energy  In  the 
force  field  are  negligible.  If,  In  the  system  under  examination,  the 
gradient  norm  Is  too  large,  the  gradient  norm  will  first  be  reduced 
using  FLEPO,  unless  LET  has  been  specified.  In  some  cases  the  FORCE 
calculation  may  be  run  only  to  decide  If  a  state  Is  a  ground  state  or  a 
transition  state.  In  which  case  the  results  have  only  two 
Interpretations.  Under  these  circumstances,  LET  may  be  warranted. 
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GRADIENT  IS  VERY  LARGE... 

In  a  calculation  of  the  thermodynamic  properties  of  the  system.  If 
the  rotation  and  translation  vibrations  are  non-zero,  as  would  be  the 
case  If  the  gradient  norm  was  significant,  then  these  "vibrations"  would 
Interfere  with  the  low-lying  genuine  vibrations.  The  criteria  for 
THERMO  are  much  more  stringent  than  for  a  vibrational  frequency 
calculation,  as  It  Is  the  lowest  few  genuine  vibrations  that  determine 
the  Internal  vibrational  energy,  entropy,  etc. 

GRADIENTS  OF  OLD  GEOMETRY,  GNORM-  nn.nnnn 

This  Is  an  error  message.  For  the  D.F.P.  geometry  optimization  to 
work  efficiently  the  gradient  norm  must  not  rise  very  much  on  making  a 
small  displacement  to  the  geometry  In  the  direction  calculated  to  reduce 
the  gradient  norm.  If  this  happens,  then  there  Is  a  finite  chance  that 
the  D.F.P.  will  fall  to  give  good  results.  Normally  the  cause  of  the 
error  Is  a  badly  specified  geometry,  and  Information  Is  printed  to  allow 
the  user  to  correct  the  fault.  If  the  user  Is  satisfied  that  there  Is 
no  fault,  then  the  key-word  "GEO-OK"  will  allow  the  calculation  to 
proceed,  but  great  care  must  then  be  taken  in  Insuring  that  the  results 
are  satisfactory.  In  particular,  verify  that  the  gradients  have  been 
reduced  to  about  zero  by  the  D.F.P.  at  the  end  of  the  calculation. 

IMPOSSIBLE  NUMBER  OF  OPEN  SHELL  ELECTRONS 

The  keyword  0PEN(nl,n2)  has  been  used,  but  for  an  even-electron 
system  nl  was  specified  as  odd  or  for  an  odd-electron  system  nl  was 
specified  as  even.  Either  way,  there  Is  a  conflict  which  the  user  must 
resolve. 

IMPOSSIBLE  OPTION  REQUESTED 

A  general  catch-all.  This  message  will  be  printed  If  two 
Incompatible  options  are  used,  such  as  both  MINDO/3  and  AMI  being 
specified.  Check  the  key-words,  and  resolve  the  conflict. 
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INTERNAL  COORDINATES  READ  IN,  AND  CALCULATION... 

If  Internal  coordinates  are  read  In,  but  the  calculation  Is  to  be 
carried  out  using  cartesian  coordinates,  then  either  all  possible 
geometric  variables  must  be  optimized,  or  none  can  be  optimized.  If 
only  some  are  marked  for  optimization,  then  ambiguity  exists.  For 
example.  If  the  bond-length  of  atom  6  Is  marked  for  optimization,  but 
the  angle  Is  not,  then  when  the  conversion  to  cartesian  coordinates 
takes  place,  the  first  coordinate  becomes  the  "X"  coordinate  and  the 
second  the  "Y"  coordinate.  These  bear  no  relationship  to  the  bond 
length  or  angle.  This  Is  a  fatal  error. 

INTERNAL  COORDINATES  READ  IN,  AND  SYMMETRY... 

If  Internal  coordinates  are  read  In,  but  the  calculation  Is  to  be 
carried  out  using  cartesian  coordinates,  then  any  symmetry  relationships 
between  the  Internal  coordinates  will  not  be  reflected  In  the  cartesian 
coordinates.  For  example.  If  the  bond-lengths  of  atoms  5  and  6  are 
equal.  It  does  not  follow  that  these  atoms  have  equal  values  for  their 
"X"  coordinates.  This  Is  a  fatal  error. 

THE  LINE  MINIMISATION  FAILED  TWICE  IN  A  ROW 

This  Is  usually  found  In  more  exotic  calculations.  For  some 
reason,  the  heat  of  formation  at  the  start  of  a  line  minimization  within 
FLEPO  was  lower  than  any  point  on  the  line  minimization,  a  situation 
that  cannot  occur  naturally. 

Look  at  the  gradient.  If  printed.  If  It  Is  acceptably  small. 
Ignore  the  message.  The  gradient  will  not  be  printed  If  It  was  very 
small:  less  than  2.0.  The  time  to  worry  Is  when  the  gradient  norm  Is 
larger  than,  say,  20:  If  larger  than  50  something  Is  badly  wrong,  and 
you  should  then  suspect  MOPAC  Is  faulty  In  some  way. 

****  MAX.  NUMBER  OF  ATOMS  ALLOWED:.... 

At  compile  time  the  maximum  sizes  of  the  arrays  In  MOPAC  are  fixed. 
The  system  being  run  exceeds  the  maximum  number  of  atoms  allowed.  To 
rectify  this,  modify  the  file  DIMSIZES.DAT  to  Increase  the  number  of 
heavy  and  light  atoms  allowed.  If  DIMSIZES.DAT  Is  altered,  then  the 
whole  of  MOPAC  should  be  re-complled  and  re-linked. 

****  MAX.  NUMBER  OF  ORBITALS:.... 

At  compile  time  the  maximum  sizes  of  the  arrays  in  MOPAC  are  fixed. 
The  system  being  run  exceeds  the  maximum  number  of  orbitals  allowed.  To 
rectify  this,  modify  the  file  DIMSIZES.DAT  to  change  the  number  of  heavy 
and  light  atoms  allowed.  If  DIMSIZES.DAT  Is  altered,  then  the  whole  of 
MOPAC  should  be  re-compiled  and  re-linked. 
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****  MAX.  NUMBER  OF  TWO  ELECTRON  INTEGRALS.. 

At  compile  time  the  maximum  sizes  of  the  arrays  In  MOPAC  are  fixed. 
The  system  being  run  exceeds  the  maximum  number  of  two-electron 
Integrals  allowed.  To  rectify  this,  modify  the  file  DIMSIZES.DAT  to 
modify  the  number  of  heavy  and  light  atoms  allowed.  If  DIMSIZES.DAT  Is 
altered,  then  the  whole  of  MOPAC  should  be  re-compiled  and  re-linked. 

NO  POINT  LOWER  IN  ENERGY  ... 

This  Is  a  non-fatal  error  message,  caused  by  the  Powell  line-search 
routine  In  the  D.F.P.  falling  to  lower  the  energy.  Two  explanations 
are  possible:  first,  the  geometry  might  be  almost  fully  optimized.  In 
which  case  the  smallest  step  possible  In  the  line  search  Is  still  too 
big;  and  second,  the  calculated  heat  of  formation  of  a  previous  step  was 
too  low.  This  would  be  an  algorithmic  error. 

This  error  message  occurs  rarely  In  routine  use,  but  can  easily  be 
Induced  If  program  parameters  such  as  the  SCF  criterion  are  altered. 

SINCE  COS  .IAS  JUST  BEEN  RESET... 

A  non-fatal  error  message.  After  two  attempts  to  lower  the  heat  of 
formation  have  failed,  the  D.F.P.  optimization  Is  terminated.  If  the 
gradients  are  acceptable,  no  further  action  Is  necessary;  otherwise 
examine  the  geometry  to  see  If  a  more  reasonable  numbering  system  can  be 
made. 


A  rarer  cause  of  this  message  Is  when  the  geometry  supplied  Is 
already  almost  fully  optimized,  and  It  becomes  optimized  within  the 
first  cycle. 

SYSTEM  DOES  NOT  APPEAR  TO  BE  OPTIMISABLE 

This  Is  a  gradient  norm  minimization  message.  These  routines  will 
only  work  If  the  nearest  minimum  to  the  supplied  geometry  In 
gradient-norm  space  Is  a  transition  state  or  a  ground  state.  Gradient 
norm  space  can  be  visualized  as  the  space  of  the  scalar  of  the 
derivative  of  the  energy  space  with  respect  to  geometry.  To  a  first 
approximation,  there  are  twice  as  many  minima  In  gradient  norm  space  as 
there  are  In  energy  space. 

It  Is  unlikely  that  there  exists  any  simple  way  to  refine  a 
geometry  that  results  In  this  message.  While  It  Is  appreciated  that  a 
large  amount  of  effort  has  probably  already  been  expended  In  getting  to 
this  point,  users  should  steel  themselves  to  writing  off  the  whole 
geometry.  It  Is  not  recommended  that  a  minor  change  be  made  to  the 
geometry  and  the  job  re-submitted. 
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TEMPERATURE  RANGE  STARTS  TOO  LOW,... 

The  thermodynamics  calculation  assumes  that  the  statistical 
summations  can  be  replaced  by  Integrals.  This  assumption  Is  only  valid 
above  100K,  so  the  lower  temperature  bound  Is  set  to  100,  and  the 
calculation  continued. 

THERE  HAVE  BEEN  3  ATTEMPTS... 

In  FLEPO,  the  energy  minimization  routine,  three  attempts  have  been 
made  to  lower  the  heat  of  formation,  and  thus  the  gradient,  but  while 
the  heat  of  formation  has  reached  a  fairly  constant  value,  the  gradient 
norm  Is  still  large.  This  Is  a  non-fatal  error,  and  almost  certainly 
the  results  will  be  useful.  Try  using  "XYZ",  or  re-arranging  the 
connectivity. 

THERE  IS  A  RISK  OF  INFINITE  LOOPING... 

The  SCF  criterion  has  been  reset  by  the  user,  and  the  new  value  Is 
so  small  that  the  SCF  test  may  never  be  satisfied.  This  Is  a  case  of 
user  beware! 

THIS  MESSAGE  SHOULD  NEVER  APPEAR,  CONSULT  A  PROGRAMMER! 

This  message  should  never  appear;  a  fault  has  been  Introduced  Into 
MOPAC,  most  probably  as  a  result  of  a  programming  error. 
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THREE  ATOMS  BEING  USED  TO  DEFINE.... 

If  the  cartesian  coordinates  of  an  atom  depend  on  the  dihedral 
angle  It  makes  with  three  other  atoms,  and  those  three  atoms  fall  In  an 
almost  straight  line,  then  a  small  change  In  the  cartesian  coordinates 
of  one  of  those  three  atoms  can  cause  a  large  change  In  Its  position. 
This  Is  a  potential  source  of  trouble,  and  the  data  should  be  changed  to 
make  the  geometric  specification  of  the  atom  In  question  less  ambiguous. 

This  message  can  appear  at  any  time,  particularly  In  reaction  path 
and  saddle-point  calculations. 

An  exception  to  this  rule  Is  If  the  three  atoms  fall  Into  an 
exactly  straight  line.  For  example,  If,  In  propyne,  the  hydrogens  are 
defined  In  terms  of  the  three  carbon  atoms,  then  no  error  will  be 
flagged.  In  such  a  system  the  three  atoms  In  the  straight  line  must  not 
have  the  angle  between  them  optimized,  as  the  finite  step  in  the 
derivative  calculation  would  displace  one  atom  off  the  straight  line  and 
the  error-trap  would  take  effect. 

Correction  Involves  re-defining  the  connectivity.  LET  and  GEO-OK 
will  not  allow  the  calculation  to  proceed. 

-  TIME  UP  - 

The  time  defined  on  the  key-words  line  or  3,600  seconds.  If  no  time 
was  specified.  Is  likely  to  be  exceeded  If  another  cycle  of  calculation 
were  to  be  performed.  A  controlled  termination  of  the  run  would  follow 
this  message.  The  job  may  terminate  earlier  than  expected:  this  Is 
ordinarily  due  to  one  of  the  recently  completed  cycles  taking  unusually 
long,  and  the  safety  margin  has  been  Increased  to  allow  for  the 
possibility  that  the  next  cycle  might  also  run  for  much  longer  than 
expected. 


TRIPLET  SPECIFIED  WITH  ODD  NUMBER  OF  ELECTRONS. 

If  TRIPLET  has  been  specified  the  number  of  electrons  must  be  even. 
Check  the  charge  on  the  system,  the  empirical  formula,  and  whether 
TRIPLET  was  Intended. 
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"""""""‘"'""'""'UNABLE  TO  ACHIEVE  SELF-CONSISTENCY 

See  the  error-message: 

«« - ****  FAILED  TO  ACHIEVE  SCF.  **** - >>» 

UNDEFINED  SYMMETRY  FUNCTION  USED 

Symmetry  operations  are  restricted  to  those  defined,  i.e.  In  the 
range  1-18.  Any  other  symmetry  operations  will  trip  this  fatal  message. 

UNRECOGNISED  ELEMENT  NAME 

In  the  geometric  specification  a  chemical  symbol  which  does  not 
correspond  to  any  known  element  has  been  used.  The  error  lies  In  the 
first  datum  on  a  line  of  geometric  data. 

****  WARNING  **** 

Don’t  pay  too  much  attention  to  this  message.  Thermodynamics 
calculations  require  a  higher  precision  than  vibrational  frequency 
calculations.  In  particular,  the  gradient  norm  should  be  very  small. 
However,  It  Is  frequently  not  practical  to  reduce  the  gradient  norm 
further,  and  to  date  no-one  has  determined  just  how  slack  the  gradient 
criterion  can  be  before  unacceptable  errors  appear  In  the  thermodynamic 
quantities.  The  0.4  gradient  norm  Is  only  a  suggestion. 

WARNING:  INTERNAL  COORDINATES. . . 

Triatomics  are,  by  definition,  defined  In  terms  of  Internal 
coordinates.  This  warning  Is  only  a  reminder.  For  diatomics,  cartesian 
and  Internal  coordinates  are  the  same.  For  tetra-atomlcs  and  higher, 
the  presence  or  absence  of  a  connectivity  table  distinguishes  Internal 
and  cartesian  coordinates,  but  for  trlatomlcs  there  Is  an  ambiguity.  To 
resolve  this,  cartesian  coordinates  are  not  allowed  for  the  data  Input 
for  trlatomlcs. 
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CRITERIA 


MOPAC  uses  various  criteria  which  control  the  precision  of  Its 
stages.  These  criteria  are  chosen  as  the  best  compromise  between  speed 
and  acceptable  errors  In  the  results.  The  user  can  override  the  default 
settings  by  use  of  key-words;  however,  care  should  be  exercised  as 
Increasing  a  criterion  can  Introduce  the  potential  for  Infinite  loops, 
and  decreasing  a  criterion  can  result  In  unacceptably  Imprecise  results. 
These  are  usually  characterized  by  'noise'  In  a  reaction  path,  or  large 
values  for  the  trivial  vibrations  In  a  force  calculation. 


9.1  SCF  CRITERION 


Name  : 

Defined  In 
Default  value 
Basic  Test 


SCFCRT. 

ITER. 

0.00001  kcal/mole 

Change  In  energy  In  kcal/mole  on  successive 
Iterations  Is  less  than  SCFCRT. 


Exceptions:  If  PRECISE  Is  specified,  SCFCRT-0. 0000001 

If  non-varlatlonal  wavefunctlon  SCFCRT-0. 0000001 
If  a  polarization  calculation 
or  gradient  minimization  SCFCRT-0. 0000001 

If  SCFCRT-n.nnn  Is  specified  SCFCRT-n.nnn 


Secondary  tests:  (1)  Change  In  density  matrix  elements  on  two 

successive  Iterations  must  be  less  than  0.001 
(2)  Change  In  energy  In  eV  on  three  successive 
Iterations  must  be  less  than  10  x  SCFCRT. 
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9.2  GEOMETRIC  OPTIMIZATION  CRITERIA 


Name : 

Defined  In 
Default  value 
8as1c  Test 

TOLERX  "Test  on  X  Satisfied" 

FLEPO 

0.0001 

The  projected  change  In  geometry  Is  less  than 

TOLERX  Angstroms. 

Exceptions 

If  PRECISE  Is  specified,  TOLERX-  0.00001 

Name: 

Defined  In 
Default  value 
Basic  Test 

EYEAD  "Herbert's  Test  Satisfied" 

FLEPO 

0.001 

The  projected  decrease  In  energy  Is  less  than 

EYEAD  Kcal s/mol e. 

Exceptions 

If  PRECISE  Is  specified,  EYEAD-0. 00001 

Name  : 

Defined  In 
Default  value 
Basic  Test 

TOLERG  "Test  on  Gradient  Satisfied" 

FLEPO 

1.0 

The  gradient  norm  In  Kcal s/mol e/Angstrom  Is 
less  than  TOLERG. 

Exceptions 

If  PRECISE  Is  specified,  T0LERG-0.002 

If  GNORM=n.nnn  Is  specified,  TOLERG-n.nnn 

If  a  SAODLE  calculation,  TOLERG  Is  made  a  function 
of  the  last  gradient  norm. 

Name: 

Defined  In 
Default  value 
Basic  Test 

TOLERF  "Heat  of  Formation  Test  Satisfied" 

FLEPO 

0.002 

The  calculated  heats  of  formation  on  two  successive 
cycles  differ  by  less  than  TOLERF. 

Exceptions 

If  PRECISE  Is  specified,  TOLERF-O. 00004 

Secondary  Tests 

For  the  TOLERG,  TOLERF,  and  TOLERX  tests,  a 
second  test  in  which  individual  component  of  the 
gradient  should  be  larger  than  TOLERG  divided  by  the 
square  root  of  the  number  of  variables  must  be 
satisfied. 

Other  Tests 

If,  after  the  TOLERG,  TOLERF,  or  TOLERX  test  has  been 
satisfied  three  consecutive  times,  and  the  heat  of 
formation  has  dropped  by  less  than  0.3Kcal/mole,  the 
optimization  Is  stopped. 

Name: 

Defined  in 
Default  Value 
Basic  Test 

T0L2 

NLLSQ 

0.000  000  000  001 

The  square  of  the  ratio  of  the  projected  change  In  the 
geometry  to  the  actual  geometry  Is  less  than  T0LS1 . 
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CHAPTER  10 


DEBUGGING 


There  are  three  potential  sources  of  difficulty  In  using  MOPAC, 
each  of  which  requires  special  attention.  There  can  be  problems  with 
data,  due  to  errors  In  the  data,  or  MOPAC  may  be  called  upon  to  do 
calculations  for  which  It  was  not  designed.  There  are  intrinsic  errors 
In  MOPAC  which  extensive  testing  has  not  yet  revealed,  but  which  a 
user's  novel  calculation  uncovers.  Finally  there  can  be  bugs  Introduced 
by  the  user  modifying  MOPAC,  either  to  make  It  compatible  with  the  host 
computer,  or  to  implement  local  features. 

For  whatever  reason,  the  user  may  need  to  have  access  to  more 
Information  than  the  normal  key-words  can  provide,  and  a  second  set, 
specifically  for  debugging.  Is  provided.  These  key-words  give 
Information  about  the  working  of  Individual  subroutines,  and  do  not 
affect  the  course  of  the  calculation. 


10.1  DEBUGGING  KEY-WORDS 

FULL  LIST  OF  KEY-WORDS  FOR  DEBUGGING  SUBROUTINES 
INFORMATION  PRINTED 


1ELEC  the  one-electron  matrix.  Note  1 

COMPFG  Heat  of  Formation. 

DCART  Cartesian  derivatives. 

DEBUG  Note  2 

DEBUGPULAY  Pulay  matrix,  vector,  and  error-function.  Note  3 

DENSITY  Every  density  matrix.  Note  1 


DERIV  All  gradients,  and  other  data  In  DERIV. 

DFORCE  Print  Force  Matrix. 

EIGS  All  eigenvalues. 

FLEPO  Details  of  D.F.P.  minimization. 

FMAT 

FOCK  Every  Fock  matrix  Note  1 

HCORE  The  one  electron  matrix,  and  two  electron  Integrals. 

ITER  Values  of  variables  and  constants  In  ITER. 

LINMIN  Function  values,  step  sizes  at  all  points  In  LINMIN. 

LOCMIN  Function  values,  step  sizes  at  a’l  points  In  LOCMIN. 

MOLDAT  Molecular  data,  number  of  orbitals,  "U"  values,  etc. 
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MECI 

C.I.  matrices,  M.O.  Indices,  etc. 

PL 

Differences  between  density  matrix  elements 

In  ITER. 

Note  4 

TIMES 

Times  of  stages  within  ITER. 

VECTORS 

All  eigenvectors  on  every  iteration. 

Note  1 

Note  the 

space  before  PL:  this  space  Is  obligatory. 

NOTES 

1.  These  key-words  are  activated  by  the  key-word  DEBUG.  Thus  If 
DEBUG  and  FOCK  are  both  specified,  every  Fock  matrix  on  every 
Iteration  will  be  printed. 

2.  DEBUG  Is  not  Intended  to  Increase  the  output,  but  does  allow 
other  key-words  to  have  a  special  meaning. 

3.  PULAY  Is  already  a  key-word,  so  DEBUGPULAY  was  an  obvious 
alternative. 

4.  PL  Initiates  the  output  of  the  value  of  the  largest  difference 
between  any  two  density  matrix  elements  on  two  consecutive 
Iterations.  This  Is  very  useful  when  Investigating  options  for 
Increasing  the  rate  of  convergance  of  the  SCF  calculation. 


SUGGESTED  PROCEDURE  FOR  LOCATING  BUGS 

Users  are  supplied  with  the  source  code  for  MOPAC,  and,  while  the 
original  code  Is  fairly  bug-free,  after  It  has  been  modified  there  Is  a 
possibility  that  bugs  may  have  been  Introduced.  In  these  circumstances 
the  author  of  the  changes  Is  obviously  responsible  for  removing  the 
offending  bug,  and  the  following  Ideas  might  prove  useful  In  this 
context. 

First  of  all,  and  most  Important,  before  any  modifications  are  done 
a  back-up  copy  of  the  standard  MOPAC  should  be  made.  This  will  prove 
Invaluable  In  pinpointing  deviations  from  the  standard  working.  This 
point  cannot  be  over-emphasized  -  MAKE  A  BACK-UP  BEFORE  MODIFYING 
MOPAC! ! ! ! 

Clearly,  a  bug  can  occur  almost  anywhere,  and  a  logical  search 
sequence  Is  necessary  In  order  to  minimize  the  time  taken  to  locate  It. 

If  possible,  perform  the  debugging  with  a  small  molecule.  In  order 
to  save  time  (debugging  Is,  of  necessity,  time  consuming)  and  to 
minimize  output. 

The  two  sets  of  subroutines  In  MOPAC,  those  Involved  with  the 
electronics  and  those  Involved  In  the  geometries,  are  kept  strictly 
separate,  so  the  first  question  to  be  answered  Is  which  set  contains  the 
bug.  If  the  heats  of  formation,  derivatives,  I.P.s,  and  charges,  etc., 
are  correct,  the  bug  lies  In  the  geometries;  If  faulty,  In  the 
electronics. 
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Bug  In  the  Electronics  Subroutines. 

Use  formaldehyde  for  this  test.  The  supplied  data-flle  MNRSD1.DAT 
could  be  used  as  a  template  for  this  operation.  Use  key-words  1SCF, 
DEBUG,  and  any  others  necessary. 

The  main  steps  are: 

(1)  Check  the  starting  one-electron  matrix  and  two-electron 
Integral  string,  using  the  key-word  HCORE.  It  Is  normally  sufficient  to 
verify  that  the  two  hydrogen  atoms  are  equivalent,  and  that  the  pi 
system  Involves  only  pz  on  oxygen  and  carbon.  Note  that  numerical 
values  are  not  checked,  but  only  relative  values. 

If  an  error  Is  found,  use  MOLDAT  to  verify  the  orbital  character, 

etc. 


If  faulty  the  error  lies  In  READ,  GETGEO  or  MOLDAT. 

Otherwise  the  error  lies  In  HCORE,  HI  ELEC  or  ROTATE. 

If  the  starting  matrices  are  correct,  go  on  to  step  (2). 

(2)  Check  the  density  or  Fock  matrix  on  every  Iteration,  with  the 
words  FOCK  or  DENSITY.  Check  the  equivalence  of  the  two  hydrogen  atoms, 
and  the  pi  system,  as  In  (1). 

If  an  error  Is  found,  check  the  first  Fock  matrix.  If  faulty,  the 
bug  lies  In  ITER,  probably  In  the  Fock  subroutines  FOCK1  or  FOCK2 .  or 
In  the  (guessed)  density  matrix  (MOLDAT).  An  exception  Is  In  the  UHF 
closed-shell  calculation,  where  a  small  asymmetry  Is  Introduced  to 
Initiate  the  separation  of  the  alpha  and  beta  UHF  wavefunctlons . 

If  no  error  Is  found,  check  the  second  Fock  matrix.  If  faulty,  the 
error  lies  In  the  density  matrix  DENSIT,  or  the  diagonal Izatlon  RSP. 

If  the  Fock  matrix  Is  acceptable,  check  all  the  Fock  matrices.  If 
the  error  starts  In  Iterations  2  to  4,  the  error  probably  lies  In  CNVG , 
If  after  that.  In  PULAY,  If  used. 

If  SCF  Is  achieved,  and  the  heat  of  formation  Is  faulty,  check 
HELECT.  If  C.I.  was  used  check  MECI. 

If  the  derivatives  are  faulty,  use  DCART  to  verify  the  cartesian 
derivatives.  If  these  are  faulty,  check  DCART  and  DHC.  If  they  are 
correct,  or  not  calculated,  check  the  DERIV  finite  difference 
calculation. 

If  the  geometric  calculation  Is  faulty,  use  FLEPO  and  LINMIN  to 
monitor  the  optimization,  DERIV  may  also  be  useful  here. 

For  the  FORCE  calculation,  DCART  or  DERIV  are  useful  for 
variational ly  optimized  functions,  COMPFG  for  non-varlatlonal ly 
optimized  functions. 
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For  reaction  paths,  verify  that  FLEPO  Is  working  correctly;  If  so, 
then  PATHS  Is  faulty. 

For  saddle-point  calculations,  verify  that  FLEPO  Is  working 
correctly;  If  so,  then  REACH  Is  faulty. 

Keep  In  mind  the  fact  that  HOPAC  Is  a  large  calculation,  and  while 
Intended  to  be  versatile,  many  combinations  of  options  have  not  been 
tested.  If  a  bug  Is  found  In  the  original  code,  please  communicate 
details  to  the  Academy,  to  Dr.  James  J.  P.  Stewart,  Frank  J.  Seller 
Research  Laboratory,  U.S.  Air  Force  Academy,  Colorado  Springs,  CO  80840 


CHAPTER  11 
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Installing  MOPAC  on  the  IBM-PC  XT  or  AT 

This  discussion  assumes  that  the  user  has  an  IBM-PC  XT  or  AT 
computer  (hereafter  referred  to  as  an  IBM-PC)  or  IBM-compatible  with  an 
8087  coprocessor,  a  formatted  10  Mb  hard  disk  (drive  C:>,  and  a  360  Kb 
floppy  disk  (drive  A:).  The  directory  structure  on  the  hard  disk  should 
contain  a  directory  for  executable  files  (where  the  FORTRAN  compiler 
resides,  the  default  Is  C:\BINV),  a  directory  for  libraries  (where  the 
FORTRAN  libraries  reside,  the  default  Is  C:\LIB\>,  and  a  directory  to 
hold  the  MOPAC  source  files,  the  default  being  C:\M0PAC\.  The  Microsoft 
FORTRAN  compiler  and  DOS  3.1  (or  higher)  must  also  be  on  the  system. 
Finally,  the  PATH  should  be  set  to  include  each  of  these  directories. 

The  steps  to  be  carried  out  to  get  a  working  copy  of  MOPAC  are  as 
follows  -  don't  worry  -  a  specific  example  of  the  commands  required  for 
these  steps  are  Included  at  the  end  of  this  chapter. 

1.  Copy  all  the  files  on  the  MOPAC  diskettes  to  the  PC  -  the  files 

were  copied  directly  onto  the  source  diskettes.  They  can  be 

retrieved  by  first  entering  the  MOPAC  directory  and  then  typing 

COPY  A:*.*  C:  for  each  diskette. 

At  this  point  you  should  have  the  following  sets  of  files: 

1.  A  set  of  73  FORTRAN  files  with  names  <f 1 lename>.F0R 

2.  A  set  of  10  dummy  FORTRAN  files  with  names  S<f11ename>.F0R 

3.  A  set  of  7  CREATE  files  with  names  <f 1 lename> .CRE  (MOPAC 

specific  postfix) 

4.  A  set  of  6  BATCH  files  with  names  <f 1 lename> .BAT 

5.  A  set  of  3  DATA  files  with  names  <f 1 lename> .DAT  (MOPAC 

specific  postfix) 

6.  A  FIL  file  called  MOPACLIN. FIL.  (MOPAC  specific  postfix) 
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2.  Compile  the  FORTRAN  source  files:  Utility  batch  files  called 
COMPILE.BAT  and  SCOMPILE.BAT  have  been  provided  for  this 
purpose. 

DIMSIZES.DAT 

All  the  array  sizes  used  In  the  FORTRAN  files  are  given  In 
DIMSIZES.DAT.  There  are  only  four  parameters  that  actually 
ever  need  to  be  changed.  If,  for  whatever  reason,  DIMSIZES.DAT 
needs  to  be  changed,  then  new  object  libraries  must  be  created 
and  linked  (COMPILE,  LIBRARY,  etc.),  as  modules  compiled  with 
different  DIMSIZES.DAT  will  be  Incompatible. 

The  parameters  within  DIMSIZES.DAT  that  the  user  can 
modify  are  MAXLIT,  MAXHEV,  MAXTIM  and  MAXDMP.  MAXLIT  Is 
assigned  a  value  equal  to  the  largest  number  of  hydrogen  atoms 
that  a  MOPAC  job  Is  expected  to  run,  MAXHEV  is  assigned  the 
corresponding  number  of  heavy  (non-hydrogen)  atoms.  MAXTIM  Is 
the  default  maximum  time  In  seconds  a  job  Is  allowed  to  run 
before  either  completion  or  a  restart  file  being  written. 
MAXDMP  Is  the  default  time  in  seconds  for  the  automatic  writing 
of  the  restart  files.  If  disk  space  Is  at  a  premium,  you  might 
want  to  set  MAXDMP  as  MAXDMP-999999. 

The  default  values  of  the  parameters  are 

MAXHEV  -  7;  Maximum  number  of  heavy  (non-hydrogen)  atoms. 

MAXLIT  =  7;  Maximum  number  of  light  (hydrogen)  atoms. 

MAXTIM  =3600;  Default  job  time  In  seconds  (see  also  "  T="). 

MAXDMP  =3600;  Default  time  interval  between  automatic  dumps 
(see  also  "DUMP-"). 

The  maximum  practical  size  for  the  arrays  is  MAXHEV=8  and 
MAXLIT=8.  In  this  case  MOPACLIB  can  only  be  linked  with 
MINLIB;  the  other  options  (FORCELIB,  etc)  produce  code  which  Is 
too  large  to  run  on  a  PC. 

One  can  actually  sometimes  run  more  heavy  or  light  atoms 
than  allowed  based  on  MAXHEV  and  MAXLIT.  For  example.  In  the 
case  where  MAXHEV  =  MAXLIT  =  7,  one  can  use  8  heavy  atoms.  If 

there  are  only  three  hydrogens  or  less  In  the  molecule,  or  one 

can  use  9  hydrogens  In  a  system.  If  there  are  5  or  less  heavy 
atoms.  The  general  rules  are: 

total  number  of  atoms  <=  MAXHEV  +  MAXLIT 
(#  heavy)*4  +  Of  light)  <*  (MAXHEV  +  MAXLIT)*3 


If  MAXHEV  and  MAXLIT  are  to  be  set  to  different  values  In 
DIMSIZES.DAT,  the  following  rules  apply: 

MAXHEV‘4  +  MAXLIT  <  (MAXHEV  +  MAXLIT)*3 
2*(MAXHEV*4  +  MAXLIT)  >  (MAXHEV  +  MAXLIT>*3 
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To  compile  all  the  FORTRAN,  run  the  file  COMPILE.  The 
command  for  this  Is 


COMPILE  (CR) 


|  A  series  of  <fi lename>.OBJ  files  will  be  created.  The  COMPILE 

j  command  takes  about  one  and  a  half  hours  to  execute  on  the  XT, 

|  one  hour  on  the  AT.  As  this  Is  the  slow  step  In  preparing 

MOPAC,  once  it  is  started  go  and  have  a  coffee.  Once  this  Is 
finished,  run  SCOMPILE  to  compile  the  dummy  routines,  this 
|  takes  about  5  minutes.  If  you  do  not  have  an  8087  coprocessor, 

S  remove  the  command  $N0FL0ATCALLS  from  the  file  METACOM.DAT 

*  before  compiling. 


i 


Create  the  link  libraries.  Batch  files  to  carry  out  this 
operation  are  MOPACLIB.CRE,  MINLIB.CRE,  etc.  Before  use,  the 
first  line  of  each  of  these  batch  files  must  be  changed  to  give 
the  path  to  the  library  directory.  As  supplied,  the  library 
directory  is  "C:VLIB\"  (line  1  in  the  <f 1 lename>.CRE  files). 
The  libraries  are  created  one  at  a  time  by  typing,  e.g., 

FLIB  0MOPACLIB.CRE 

(Assuming  FLIB  Is  the  name  of  your  MSFORTRAN  library  manager.) 
The  MS  FORTRAN  linker  comes  with  the  name  LINK.EXE,  as  does  the 
MSDOS  linker,  so  rename  the  FORTRAN  linker  FLINK.EXE.  The 
batch  file  assumes  libraries  with  these  names  do  not  already 
exist.  If  they  do  exist,  then  delete  them  first.  This  step 
takes  very  little  time. 


The  full  list  of  <f11ename>.CRE  files  Is  as  follows: 


(a)  MOPACLIB.CRE 

(b)  MINLIB.CRE 

(c)  ENPARTLI.CRE 

(d)  FORCELIB.CRE 

(e)  MECILIB.CRE 

(f)  NLLSQLIB.CRE 

(g)  SADDLELI.CRE 


The  main  portion  of  MOPAC. 

Linked  to  create  the  minimum  configuration. 

Linked  to  create  the  version  with  energy  partitioning 

.  FORCE  calculations. 

"  . . MECI  calculations. 

.  NLLSQ  calculations. 

"  . .  SADDLE  calculations 


If  you  want  to  create  all  the  libraries  at  once  you  can 
use  the  batch-file  LIBRARY.BAT  to  create  them.  LIBRARY  will 
also  automatically  delete  old  libraries  with  the  same  name. 
LIBRARY  assumes  these  files  are  in  the  c:\LIB\  directory.  If 
not,  modify  LIBRARY.BAT. 


4.  Create  an  executable  program.  The  linking  of  the  object  files 
and  the  creation  of  an  executable  program  are  carried  out  by 
the  batch-file  DOFLNK  and  the  command  file  MOPACLIN.FIL. 
DOFLNK  assumes  the  executable  of  your  FORTRAN  linker  Is  called 
FLINK.EXE  (The  MS  FORTRAN  linker  comes  called  LINK.EXE,  as  does 
the  MSDOS  linker.)  The  second  line  of  the  file  MOPACLIN.FIL 
gives  the  title  that  the  executable  file  will  have.  Line  three 
Indicates  the  name  given  to  the  linker  map  produced  during 
linking.  The  fourth  line  Indicates  what  libraries  will  be  used 
In  the  link  step.  By  default,  the  file  supplied  will  set  up 
the  minimum  MOPAC  configuration.  Change  "MINLIB"  to  whatever 
version  you  want  to  create  (MOPACLIB,  FORTRAN,  and  8087  are 
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always  required  -  If  you  don't  have  an  8087  processor  chip 
change  "8087"  to  "MATH"  and  remove  the  command  SNOFLOATCALLS 
from  the  file  METAC0M.DAT).  Linking  takes  about  10  minutes. 

5.  Cleaning  up.  Having  created  an  executable  version  of  MOPAC, 
most  of  the  <f 1 lena1e> .OBJ  files  can  now  be  deleted.  The 
source  files  <f1 lename>.F0R  can  be  deleted  as  well.  If  you  need 
to  save  space  on  your  disk.  The  batch  file  CLEAN.BAT  will 
delete  all  the  object  files  now  contained  In  MOPACLIB.LIB,  but 
leaving  the  ones  required  for  creating  FORCELIB.LIB,  etc.  The 
batch  file  SCLEAN.BAT  will  delete  the  remaining  <f 1 lename>.0BJ 
files  except  for  MND0.0BJ  and  BLOCK. OBJ,  which  are  needed  for 
linking  the  different  forms  of  MOPAC.  LCLEAN.BAT  will  delete 
the  list  files  created  during  library  creation  and  the  loader 
map  created  during  linking.  Either  the  executable  can  be 

transferred  to  a  directory  (\BIN\,  \EXE\,  SYS\,  etc)  which 
holds  all  your  executable  code,  or  left  where  It  Is. 


Sample  Session  for  Installing  MOPAC  on  an  IBM-PC  XT 


This  example  assumes  that  the  user  has  an  IBM  XT  or  IBM-compatible 
with  an  8087  coprocessor,  a  hard  disk  drive  (drive  C:>,  a  floppy  disk 
drive  (drive  A:),  640  K  RAM,  and  a  copy  of  MS  DOS  version  3.10  in  a 
subdirectory  on  the  hard  disk  called  C:\SYS\.  The  user  Is  also  assumed 
to  have  a  copy  of  Microsoft  FORTRAN  3.31  (other  versions  should  work), 
and  a  copy  of  the  FORTRAN  source-code  for  MOPAC  3.10  on  floppy  disks, 
both  of  which  must  be  Installed. 

On  power  up,  the  machine  should  be  In  the  root  directory  C:\.  From 
here  you  have  to  create  a  directory  for  the  FORTRAN  compiler  (C:\BIN\), 
a  directory  for  the  libraries  (C:\LIB\),  and  a  directory 
for  MOPAC  (C:\M0PACV).  Type: 

MD  BIN  (If  this  directory  doesn’t  already  exist) 

MD  LIB  (If  this  directory  doesn't  already  exist) 

MD  MOPAC 

Next  copy  the  FORTRAN  files  Into  the  proper  directories.  Insert  the 
FORTRAN  diskettes  into  drive  A:  and  type: 

COPY  A:*. EXE  C:\BIN 
COPY  A:*. LIB  C:\LIB 
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Next  change  the  names  of  the  FORTRAN  linker  and  library  programs. 

RENAME  C:\BIN\LINK.EXE  C:\BIN\FLINK.EXE 
RENAME  C:\BIN\LIB.EXE  C:\BIN\FLIB.EXE 

This  concludes  putting  FORTRAN  on  the  machine.  If  FORTRAN  was  already 
on  the  machine,  the  only  directory  that  had  to  be  made  was  MOPAC.  The 

default  labels  In  the  .BAT  files  must  be  changed  to  Indicate  where  your 

libraries  reside  If  the  are  not  In  C:\LIB\  and  the  path  must  be  set 

(using  SET  PATH)  to  your  FORTRAN  executable  files. 

The  next  step  Is  to  create  (or  alter)  the  file  AUTOEXEC.BAT  In  the 
root  directory.  This  session  will  assume  you  have  to  use  the  EDLIN 
editor  (part  of  MS  DOS)  to  edit  and  that  AUTOEXEC.BAT  doesn't  exist. 

To  create  It  type: 

EDLIN  AUTOEXEC.BAT 
I 

SET  PATH-C:\;C:\SYS\;C:\BIN\ 

SET  LIB-C: \LIB\ 

<ctrlxbreak> 

E 

If  you  keep  your  MS  DOS  files  in  a  directory  other  than  C:\SYS\,  modify 
the  SET  PATH  to  reflect  this.  If  you  keep  your  FORTRAN  somewhere  other 
than  C:\BINN,  again  modify  the  PATH.  If  your  libraries  are  kept  In 
other  than  the  C:\LIB\  subdirectory,  modify  the  SET  LIB  command.  If 
AUTOEXEC.BAT  exists,  use  EDLIN  to  add  or  modify  these  commands. 

Now  you  are  ready  to  Install  MOPAC.  Type: 

AUTOEXEC 

to  run  your  new  AUTOEXEC  file,  then  enter  the  MOPAC  subdirectory: 

CD  \ MOPAC 

For  each  of  the  MOPAC  diskettes  -  put  them  In  drive  A:  and  type: 

COPY  A:*.*  C: 

Then  produce  the  executable. 

COMPILE 

SCOMPILE 

LIBRARY 

DOFLNK 

These  commands  will  produce  the  minimum  MOPAC  configuration  (no  FORCE, 
NLLSQ,  etc.).  To  produce  a  version  with  one  of  these  keywords  (only  one 
Is  possible!)  replace  the  reference  to  MINLIB  In  MOPACLIN.FIL  with  the 
appropriate  library  (FORCELIB,  NLLSQLIB,  etc.).  You  will  also  need  to 
change  the  two  references  to  the  name  of  the  executable  as  well.  For 
example  change  MOPACMIN  to  MOPACFOR  when  creating  an  executable  capable 
of  running  FORCE  calculations. 
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Suggested  names  of  executables  are: 


MOPACLIB  linked  with 

MINLIB 

ENPARTLI 

FORCELIB 

MECILIB 

NLLSQLIB 

SADDLELI 


Names  of  Executable 

MOPACMIN 

MOPACENP 

MOPACFOR 

MOPACMEC 

MOPACNLL 

MOPACSAD 


You  don't  have  to  use  these  names  for  the  executables,  but  they  are 
recommended. 

To  clean  off  files  no  longer  needed  type: 

CLEAN  (object  files  used  to  create  MOPACLIB  deleted) 

SCLEAN  (object  files  used  to  create  MINLIB,  FORCELIB,  etc. 
deleted) 

This  concludes  the  sample  session. 


How  to  use  MOPAC 


A  batch  file  MOPAC.BAT  Is  provided  to  run  MOPAC.  It  assumes  that 
an  executable  called  <executab1e>.exe  exists  and  that  a  data-flle  called 
<f1 lename>.DAT  exists. 

The  command  to  run  MOPAC  is  as  follows: 

MOPAC  <f11ename>  <executable> 
or 

MOPAC  <f11ename> 
e.g. 

MOPAC  MNRSD1  MOPACFOR 

If  the  name  of  the  executable  Is  left  out,  the  batch  file  defaults 
to  the  executable  MOPACMIN.EXE. 

While  the  job  Is  running.  Information  about  each  SCF  Iteration  and 
times  for  each  SCF  cycle  will  be  displayed  on  the  screen;  the  rest  of 
the  output  Is  directed  to  files.  If  the  job  crashes  (by  typing  (Cntl)-C 
(Ret)  Y  (Ret))  after  a  DUMP  and  before  completion,  the  following  files 
will  exist: 

F0R009  Restart  file 
F0R010  Density  file 

These  need  to  be  renamed  to  allow  the  RESTART  to  work.  Rename  them  as 
fol lows: 

Rename  "F0R009"  "<f 1 lename>.RES" . 
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Rename  "FOROIO"  "<f11ename>.DEN". 

If  the  job  stops  before  completion  due  to  running  out  of  time,  the 
restart  and  output  files 

<f11ename>.RES  Restart  file 
<f1 lename>.OEN  Density  file 
<f1 lename>.OUT  Output  file 

will  be  created.  Finally,  If  the  job  runs  to  completion,  the  archive 
file 

<f1 lename>.ARC  Archive  file 
will  also  be  made. 

For  long  jobs  the  facility  exists  to  stop  the  job  while  It  Is 
running.  Every  so  often,  about  once  every  60  minutes,  a  restart  file  Is 
written  and  closed.  If  you  run  a  job  overnight,  say,  and  want  to  use 
the  PC  for  other  work,  then  crash  the  current  job  (  (Cntl)-C  (Ret)  then 
answer  Y  (Ret)  to  the  prompt  -  do  you  want  to  terminate  the  batch  job?). 
Rename  F0R009  and  FOROIO  as  described  above,  then  modify  the  data-flle 
to  Include  "RESTART",  and  resubmit  It  when  ready.  It  will  pick  up  where 
it  left  off,  with  an  average  loss  of  time  of  30  minutes.  Of  course,  the 
same  strategy  will  work  If  there  Is  a  power  failure. 

With  the  keyword  RESTART,  restarts  should  be  transparent.  When 
MOPAC  creates  restart  files,  do  not  change  them  -  they  are  required  "as 
Is"  for  the  job  to  restart. 
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Size  of  MOPAC 

The  following  sizes  of  MOPAC  modules  Is  provided  to  give  users  a 
feel  for  how  big  various  parts  are. 

MOPACMIN  Code  only  228  Kb 

MOPACENP  244  Kb 

MOPACMEC  234  Kb 

MOPACNLL  220  Kb 

MOPACFOR  227  Kb 

MOPACSAD  237  Kb 

MOPACMIN  with  1  heavy  and  1  light  atom  379  Kb 

MOPACMIN  with  7  heavy  and  7  light  atoms  536  Kb 

MOPACMIN  with  8  heavy  and  8  light  atoms  588  Kb 

MOPACMIN  with  9  heavy  and  9  light  atoms  640  Kb 

Approximate  size  of  MOPACMIN  with  n  heavy  and  n  light  atoms  is 
376  +  3.26  *  n  *  n  Kb 
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AABABC 

CHRGE 

DENSIT 

DIAT 

FLEPO 

FORCE 

GETGEO 

HCORE 

LINMIN 

MECI 

NLLSQ 

PERM 

REFER 

SPCG 

WRITE 


FORTRAN  FILES 


TABLE  1 

NAMES  OF 

FORTRAN-77 

ANAVIB 

AXIS 

CNVG 

COMPFG 

DEPVAR 

DERIV 

DIPOLE 

DOT 

FMAT 

FOCK1 

FORSAV 

FRAME 

GETSYM 

GMETRY 

HELECT 

IJKL 

LOCAL 

LOCMIN 

MNDO 

MOLDAT 

NUCHAR 

OSINV 

PULAY 

REACT 1 

REPP 

ROTATE 

SWAP 

SYMTRY 

WRTKEY 

XYZINT 

FILES 
BLOCK 
DCART 
DFPSAV 
ENPART 
FOCK2 
FREQCY 
HI  ELEC 
INDEX 
MAMULT 
MULLIK 
PARSAV 
READ 
RSP 

THERMO 
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SUBROUTINE  CALLS  IN  MOPAC 

A  list  of  the  program  segments  which  call  various  subroutines. 
Note  that  function  calls  are  not  Included  here! 

Note  Subroutine  Called  by 

AABABC  MECI 

AABACD  MECI 

AABBCD  MECI 

ANAVIB  FORCE 

AXIS  FORCE  FRAME 

BABBBC  MECI 

BABBCD  MECI 

BANGLE  XYZINT 

BFN  DIAT 

a  BINTGS  DIAT2 

BONDS  WRITE 

a  CALPAR  AMI 

CHRGE  FMAT  FOCK1  WRITE 

a  CNVG  ITER 

COMPFG  FLEPO  FMAT  FORCE  LINMIN  MNDO 

REACT 1  LOCMIN  NLLSQ 
a  COE  DIAT 

a  DANG  XYZINT 

a  DATE  WRITE 
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DCART 

DERIV 

a 

DENROT 

WRITE 

DENSIT 

ITER  MULLIK 

DEPVAR 

HADDON 

DERIV 

COMPFG  WRITE 

a 

DFPSAV 

FLEPO  PATHS 

a 

DHC 

DERIV 

DIAG 

ITER 

a 

DIAGI 

MECI 

DIAT 

HI  ELEC 

a 

DIHED 

XYZINT 

DIPOLE 

WRITE  FMAT 

DOT 

FLEPO  FORCE  LOCMIN  NLLSQ 
SEARCH  WRITE 

PULAY 

REACT 1 

a 

ENPART 

WRITE 

EXCHNG 

LINMIN  LOCMIN 

r 

b 

FLEPO 

MNDO  PATHS  REACT 1  FORCE 

FMAT 

FORCE 

b 

FOCK1 

ITER 

b 

FOCK2 

DCART  ITER 

b 

FOCK2D 

DCART 

FORCE 

MNDO 

a 

FORSAV 

FMAT 

a 

FRAME 

FMAT  FORCE  FREQCY 

FREQCY 

FORCE 

GEOUT 

FLEPO  GMETRY  REACT1  READ 
NLLSQ 

WRITE 

GETGEO 

REACT1  READ 

GETSYM 


READ 


r/UMiU  I  ykMMT,  tUMMOfl  &&S&8i2s  rV»MCS _ fiGMKft  _ »ayi>:l 
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GMETRY 

COMPFG  DERIV  FMAT  FORCE  MOLDAT 
WRITE  DENROT  MULLIK  POLAR  REACT1 

GRID 

MAIN 

b 

H1ELEC 

DCART  DIAT  HCORE 

a 

HADOON 

DEPVAR  SYMTRY 

b 

HCORE 

COMPFG  DERIV 

HELECT 

DCART  ITER 

IJKL 

MECI 

b 

ITER 

COMPFG  DERIV 

LINMIN 

FLEPO 

LOCAL 

WRITE 

LOCMIN 

NLLSO 

a 

MAMULT 

PULAY 

MATOUT 

FORCE  WRITE  FMAT  ITER  LOCAL  MECI 

MECI 

ITER  WRITE 

b 

MOLDAT 

MNDO 

MULLIK 

WRITE 

MULT 

MULLIK 

NLLSO 

FORCE  MNDO 

NUCHAR 

GETSYM  READ 

OSINV 

PULAY 

a 

PARSAV 

NLLSO 

PATHS 

MNDO 

PERM 

MECI 

PULAY 

ITER 

REACT1 

MNDO 

b 

READ 

MNDO 

READA 

DEPVAR  DERIV  FLEPO  FMAT  FORCE  GETGEO 
ITER  MECI  MOLDAT  NLLSQ  NUCHAR 
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I 

! 

*’ 

J 

* 


REACT1  THERMO  WRITE 

REFER 

MOLDAT 

a 

REPP 

ROTATE 

b 

ROTATE 

DERIV  HCORE 

RSP 

FMAT  FORCE  ITER  MULLIK 

FREQCY  MECI  AXIS 

SECOND 

DERIV  FLEPO  FMAT  FORCE  ITER  MNDO 

NLLSQ  PATHS  REACT1  WRITE 

SPCG 

IJKL 

SS 

DIAT 

a 

SWAP 

ITER 

SYMTRY 

COMPFG  DERIV  READ  WRITE 

THERMO 

FORCE 

VECPRT 

BONDS  FORCE  HCORE  ITER  MECI  MOLDAT  MULLIK 
WRITE 

b 

WRITE 

FORCE  ITER  MNDO  PATHS  REACH 

a 

WRTKEY 

READ 

XYZINT 

GEOUT  DFPSAV  WRITE  FORCE 

NOTES: 

XYZGEO 

XYZINT 

a  It  Is  not  Important  to  know  how  these  subroutines  work  In  order  to 
understand  the  working  of  the  program. 


b  These  subroutines  form  the  backbone  of  the  program,  and  should  be 
among  the  first  to  be  looked  at  when  learning  how  It  works. 
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A  list  of  subroutines  called  by  various  segments  (the  Inverse  of  the 
first  list) 

SUBROUTINE  CALLS 


AABABC 

AABACD 

AABBCD 

AINTGS 

ANAVIB 

AXIS  RSP 

BONDS  VECPRT 

COMPF  COMPFG 


COMPFG 

SYMTRY 

GMETRY 

HCORE 

ITER 

DERIV 

DCART 

DHC 

DENROT 

GMETRY 

COE 

DEPVAR 

DERIV 

SYMTRY 

GMETRY 

HCORE 

ITER 

DCART 

DFPSAV 

XY2INT 

DIAT 

COE  DIAT2  BFN 

EXCHNG 

FLEPO 

DFPSAV 

COMPFG 

GEOUT 

LINMIN 

FMAT 

GMETRY 

RSP 

FORSAV 

COMPFG 

CHRGE 

COMPFG  VECPRT  FRAME 

F0CK1 

CHRGE 

FORCE 

GMETRY 

MATOUT 

COMPFG 

FREQCY 

NLLSQ 

THERMO 

WRITE 

FRAME 

FMAT  RSP 

XYZINT  ANAVIB 

FRAME 

AXIS 

FREQCY 

RSP  FRAME 

GEOUT 

XYZINT 

GETGEO 

XYZINT 
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GETSYM 

NUCHAR 

GMETRY 

GEOUT 

HI  ELEC 

DIAT 

HADDON 

DEPVAR 

HCORE 

HI  ELEC  ROTATE  VECPRT 

ITER 

VECPRT  FOCK2  F0CK1  PULAY  RSP  DIAG 
SWAP  WRITE  MATOUT 

DENSIT  CNVG 

LINMIN 

COMPFG  EXCHNG 

LOCAL 

MATOUT 

LOCMIN 

COMPF  EXCHNG 

MECI 

RSP  IJKL  PERM  VECPRT  MATOUT 

MNDO 

READ  MOLDAT  REACT1  WRITE  PATHS  NLLSQ  COMPFG  FORCE 

FLEPO 

MOLDAT 

GMETRY  VECPRT  REFER 

MULLIK 

RSP  GMETRY  MULT  DENSIT  VECPRT 

NLLSQ 

PARSAV  COMPFG  GEOUT  LOCMIN 

PATHS 

FLEPO  WRITE  DFPSAV 

PULAY 

MAMULT  OSINV 

REACT 1 

GETGEO  GEOUT  FLEPO  COMPFG  WRITE 

READ 

GETGEO  WRTKEY  GETSYM  NUCHAR  SYMTRY 

GEOUT  GMETRY 

ROTATE 

REPP 

RSP 

EPSETA  TRED2  TQLRAT  TQL2  TRBAK3 

SEARCH 

COMPFG 

SYMTRY 

VALUE 

VECPRT 

WRITE 


XYZINT 


HADDON 

COMPFG 


DATE  DERIV  VECPRT  MATOUT  GMETRY  GEOUT  CHRGE  BONDS 
LOCAL  ENPART  XYZINT  MULLIK  SYMTRY 


XYZGEO  BANGLE  DIHED  DANG 
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DESCRIPTION  OF  SUBROUTINES  IN  MOPAC 

AABABC  Utility:  Calculates  the  configuration  interaction  matrix 
element  between  two  configurations  differing  by  exactly 
one  alpha  M.O.  Called  by  MECI  only. 

AABACD  Utility:  Calculates  the  configuration  interaction  matrix 
element  between  two  configurations  differing  by  exactly 
two  alpha  M.O.s.  Called  by  MECI  only. 

AABBCD  Utility:  Calculates  the  configuration  interaction  matrix 
element  between  two  configurations  differing  by  exactly 
two  M.O.s;  one  configuration  has  alpha  M.O.  "A"  and  beta 
M.O.  "C"  while  the  other  configuration  has  alpha  M.O.  "B" 
and  beta  M.O.  "D".  Called  by  MECI  only. 

ANAVIB  Utility:  Gives  a  brief  interpretation  of  the  modes  of 
vibration  of  the  molecule.  The  principal  pairs  of  atoms 
Involved  in  each  vibration  are  identified,  and  the  mode 
of  motion  (tangential  or  radial)  is  output 

AXIS  Utility:  Works  out  the  three  principal  moments  of  inertia 
of  a  molecule.  If  the  system  Is  linear,  one  moment  of  inertia 
Is  zero.  Prints  moments  in  units  of  cm**(-l>  and 
10**<-40)  gram-cm-cm. 

BABBBC  Utility:  Calculates  the  configuration  Interaction  matrix 
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element  between  two  configurations  differing  by  exactly 
one  beta  M.O.  Called  by  HECI  only. 

BABBCD  Utility:  Calculates  the  configuration  Interaction  matrix 
element  between  two  configurations  differing  by  exactly 
two  beta  M.O.s.  Called  by  MECI  only. 

BONDS  Utility:  Evaluates  and  prints  the  valencies  of  atoms  and 
bond-orders  between  atoms.  Main  argument:  density  matrix. 

No  results  are  passed  to  the  calculation,  and  no  data 
are  changed.  Called  by  WRITE  only. 

CHRGE  Utility:  Calculates  the  total  number  of  valence  electrons 
on  each  atom.  Main  arguments:  density  matrix,  array  of 
atom  charges  (empty  on  Input).  Called  by  ITER  only. 

CNVG  Utility:  Used  In  SCF  cycle.  CNVG  does  a  three-point 
interpolation  of  the  last  three  density  matrices. 

Arguments:  Last  three  density  matrices.  Number  of  iterations, 
measure  of  self-consistency  (empty  on  Input).  Called  by  ITER 
only. 

COMPFG  Main  Sequence:  Evaluates  the  total  heat  of  formation  of  the 
supplied  geometry,  and  the  derivatives.  If  requested.  This 
Is  the  nodal  point  connecting  the  electronic  and  geometric 
parts  of  the  program.  Main  arguments:  on  Input:  geometry, 
on  output:  heat  of  formation,  gradients. 

DENSIT  Utility:  Constructs  the  Coulson  electron  density  matrix  from 
the  eigenvectors.  Main  arguments:  Eigenvectors,  No.  of  singly 


-  106  - 


✓ 


-- '  "  ■  '  ■>>>• 


**-“.••**  *  .  ■<  *  -  -  • 


DESCRIPTION  of  subroutines  in  mopac 


Page  C-3 


and  doubly  occupied  levels,  density  matrix  (empty  on  Input) 

Called  by  ITER. 

OENROT  Utility:  Converts  the  ordinary  density  matrix  Into 

a  condensed  density  matrix  over  basis  functions  s  (sigma), 
p  (sigma)  and  p  (pi),  l.e.,  three  basis  functions.  Useful 
in  hybridization  studies.  Has  capability  to  handling  "d" 
functions,  If  present. 

DEPVAR  Utility:  A  symmetry-defined  "bond  length"  Is  related  to 
another  bond  length  by  a  multiple.  This  special  symmetry 
funtlon  is  Intended  for  use  in  Cluster  calculations. 

Called  by  HADDON. 

DERIV  Main  Sequence:  Calculates  the  derivatives  of  the  energy  with 
respect  to  the  geometric  variables.  This  Is  done  either  by  using 
Initially  cartesian  derivatives  (normal  mode)  or  by  full  SCF 
calculations  (half-electron  and  C.I.  mode).  Arguments:  on  Input: 
geometry,  on  output:  derivatives.  Called  by  COMPFG. 

DFPSAV  Utility:  Saves  and  restores  data  used  by  the 

Davldon-Fletcher-Powel 1  geometry  optimization.  Main  arguments: 
parameters  being  optimized,  gradients  of  parameters,  last  heat 
of  formation.  Integer  and  real  control  data.  Called  by  FLEPO. 

DIAG  Utility:  Rapid  pseudo-diagonal Izatlon.  Given  a  set  of  vectors 
which  almost  block-diagonal Ize  a  secular  determinant,  DIAG 
modifies  the  vectors  so  that  the  block-diagonal Izatlon  Is  more 
exact.  Main  arguments:  Old  vectors,  Secular  Determinant, 

New  vectors  (on  output).  Called  by  ITER. 

DIAGI  Utility:  Calculates  the  electronic  energy  arising  from 
a  given  configuration.  Called  by  MECI. 

DIAT  Utility:  Calculates  overlap  Integrals  between  two  atoms  In 
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general  cartesian  space.  Principal  quantum  numbers  up  to  6,  and 
angular  quantum  numbers  up  to  2  are  allowed.  Main  arguments: 
Atomic  numbers  and  cartesian  coordinates  In  Angstroms  of  the 
two  atoms.  Diatomic  overlaps  (on  exit).  Called  by  H1ELEC. 

DIPOLE  Utility:  Evaluates  and.  If  requested,  prints  dipole  components 
and  dipole  for  the  molecule  or  Ion.  Arguments:  Density  matrix. 
Charges  on  every  atom,  coordinates,  dipoles  (on  exit). 

Called  by  WRITE  and  FMAT. 

DOT  Utility:  Given  two  vectors,  X  and  Y,  of  length  N,  function  DOT 
returns  with  the  dot  product  X.Y.  I.e.,  If  X-Y,  then  DOT  -  the 
square  of  X.  Called  by  FLEPO. 

ENPART  Utility:  Partitions  the  energy  of  a  molecule  Into  Its  monatomic 
and  diatomic  components.  Called  by  WRITE  when  the  keyword 
ENPART  Is  specified.  No  data  are  changed  by  this  call. 

EXCHNG  Utility:  Dedicated  procedure  for  storing  3  parameters  and  one 
array  In  a  store.  Used  by  LINMIN. 

FLEPO  Main  Sequence:  Optimizes  a  geometry  by  minimizing  the  energy. 
Makes  use  of  the  first  and  estimated  second  derivatives  to 
achieve  this  end.  Arguments:  Parameters  to  be  optimized, 
(overwritten  on  exit  with  the  optimized  parameters),  Number  of 
parameters,  final  optimized  heat  of  formation.  Called  by  MAIN, 
REACT 1 ,  and  FORCE. 

FMAT  Main  sequence:  Calculates  the  exact  Hessian  matrix  for  a  system 
This  Is  done  by  either  using  differences  of  first  derivatives 
(normal  mode)  or  by  four  full  SCF  calculations  (half  electron 
or  C.I.  mode).  Called  by  FORCE. 

FOCK1  Utility:  Adds  on  to  Fock  matrix  the  one-center  two  electron 
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terms.  Called  by  ITER  only. 

FO CK2  Utility:  Adds  on  to  Fock  matrix  the  two-center  two  electron 

terms.  Called  by  ITER  and  DERIV.  In  ITER  the  entire  Fock  matrix 
Is  filled;  In  OERIV,  only  diatomic  Fock  matrices  are 
constructed. 

F0CK2D  Utility:  Virtually  identical  to  F0CK2,  but  with  the  two-electron 
matrix  In  double  precision.  Called  by  DHC. 

FORCE  Main  sequence:  Performs  a  force-constant  and  vibrational 
frequency  calculation  on  a  given  system.  If  the  starting 
gradients  are  large,  the  geometry  Is  optimized  to  reduce  the 
gradient  norm,  unless  LET  Is  specified  In  the  key-words. 

Isotopic  substitution  Is  allowed.  Thermochemical  quantities 
are  calculated.  Called  by  MAIN. 

FORSAV  Utility:  Saves  and  restores  data  used  In  FMAT  In  FORCE 
calculation.  Called  by  FMAT. 

FRAME  Utility:  Applies  a  very  rigid  constraint  on  the  translations 
and  rotations  of  the  system.  Used  to  separate  the  trivial 
vibrations  In  a  FORCE  calculation. 

FREQCY  Main  sequence:  Final  stage  of  a  FORCE  calculation.  Evaluates 
and  prints  the  vibrational  frequencies  and  modes. 

GEOUT  Utility:  Prints  out  the  current  geometry.  Can  be  called  at 
any  time.  Does  not  change  any  data. 

GETGEO  Utility:  Reads  in  geometry  In  character  mode  from  specified 
channel,  and  stores  parameters  In  arrays.  Some  error-checking 
Is  done.  Called  by  READ  and  REACT1. 

GETSYM  Utility:  Reads  In  symmetry  data.  Used  by  READ. 

GMETRY  Utility:  Fills  the  cartesian  coordinates  array.  Data  are 
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supplied  from  the  array  GEO,  GEO  can  be  (a)  In  Internal 
coordinates,  or  <b)  In  cartesian  coordinates.  If  STEP  Is 
non-zero,  then  the  coordinates  are  modified  In  light  of  the 
other  geometry  and  STEP.  Called  by  HCORE,  DERIV,  READ,  WRITE, 
MOLDAT,  etc. 

HIELEC  Utility:  Given  any  two  atoms  in  cartesian  space,  H1ELEC 
calculates  the  one-electron  energies  of  the  off-diagonal 
elements  of  the  atomic  orbital  matrix. 

H(1,J>  -  -S(1 ,j)*(beta(1 )+beta(j>)/2. 

Called  by  HCORE  and  DERIV. 

HADDON  Utility:  The  symmetry  operation  subroutine,  HADDON  relates  two 
geometric  variables  by  making  one  a  dependent  function  of  the 
other.  Called  by  SYMTRY  only. 

HCORE  Main  sequence:  Sets  up  the  energy  terms  used  In  calculating  the 
SCF  heat  of  formation.  Calculates  the  one  and  two  electron 
matrices,  and  the  nuclear  energy. 

Called  by  COMPFG. 

HELECT  Utility:  Given  the  density  matrix,  and  the  one  electron  and 
Fock  matrices,  calculates  the  electronic  energy.  No  data  are 
changed  by  a  call  of  HELECT.  Called  by  ITER  and  DERIV. 

IJKL  Utility:  Fills  the  large  two-electron  array  over  a  M.O. 
basis  set.  Calls  SPCG,  and  Is  called  by  MECI. 

ITER  Main  sequence:  Given  the  one  and  two  electron  matrices,  ITER 
calculates  the  Fock  and  density  matrices,  and  the  electronic 
energy.  Called  by  COMPFG. 

LINMIN  Main  sequence:  Dedicated  to  FLEPO,  LINMIN  does  a  line-search 
to  minimize  the  heat  of  formation  of  a  system.  Arguments: 
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starting  geometry,  search  direction.  Initial  step  size. 

Initial  heat  of  formation;  on  exit:  optimized  geometry, 
optimized  heat  of  formation,  flags.  Called  by  FLEPO. 

LOCAL  Utility:  Given  a  set  of  occupied  eigenvectors,  produces  a 
canonical  set  of  localized  bonding  orbitals,  by  a  series  of 
2x2  rotations  which  maximize  <ps1**4>.  Called 
by  WRITE. 

LOCMIN  Main  sequence:  In  a  gradient  minimization,  LOCMIN  does  a  line- 
search  to  find  the  gradient  norm  minimum.  Main  arguments: 
current  geometry,  search  direction,  step,  current  gradient 
norm;  on  exit:  optimized  geometry,  gradient  norm. 

MAMULT  Utility:  Matrix  multiplication.  Two  matrices,  stored  as  lower 
half  triangular  packed  arrays,  are  multiplied  together,  and  the 
result  stored  in  a  third  array  as  the  lower  half  triangular 
array.  Called  from  PULAY. 

MATOUT  Utility:  Matrix  printer.  Prints  a  square  matrix,  and  a 

row-vector,  usually  eigenvectors  and  eigenvalues.  The  Indices 
printed  depend  on  the  size  of  the  matrix:  they  can  be  either 
over  orbitals,  atoms,  or  simply  numbers,  thus  M.O.'s  are  over 
orbitals,  vibrational  modes  are  over  numbers.  Called  by  WRITE, 
FORCE. 

MECI  Main  sequence:  Main  function  for  Configuration  Interaction, 

MECI  constructs  the  appropriate  C.I.  matrix,  and  evaluates  the 
roots,  which  correspond  to  the  electronic  energy  of  the  states 
of  the  system.  The  appropriate  root  Is  then  returned. 

Called  by  ITER  only. 

MNDO  Main  sequence:  MAIN  program.  MNDO  first  reads  In  data  using 


-  Ill 


DESCRIPTION  OF  SUBROUTINES  IN  MOPAC 


Page  C-8 


READ,  then  calls  either  FLEPO  to  do  geometry  optimization, 

FORCE  to  do  a  FORCE  calculation,  PATHS  for  a  reaction  with  a 
supplied  coordinate,  NLLSQ  for  a  gradient  minimization  or 
REACT1  for  locating  the  transition  state.  Starts  the  timer. 

MOLDAT  Main  Sequence:  Sets  up  all  the  Invariant  parameters  used  during 
the  calculation,  e.g.  number  of  electrons.  Initial  atomic 
orbital  populations,  number  of  open  shells,  etc.  Called  once  by 
MNDO  only. 

MULLIK  Utility:  Constructs  and  prints  the  Mulllken  Population 
Analysis.  Available  only  for  RHF  calculations.  Called  by 
WRITE. 

MULT  Utility:  Used  by  MULLIK  only,  MULT  multiplies  two  square 
matrices  together. 

NLLSQ  Main  sequence:  Used  In  the  gradient  norm  minimization. 

OSINV  Utility:  Inverts  a  square  matrix.  Called  by  PULAY  only. 

PARSAV  Utility:  Stores  and  restores  data  used  In  the  gradient-norm 
minimization  calculation. 

PATHS  Main  sequence:  Given  ^  reaction  coordinate  as  a  row-vector, 

PATHS  performs  a  FLEPO  geometry  optimization  for  each  point, 
the  later  geometries  being  Initially  guessed  from  a  knowledge 
of  the  already  optimized  geometries,  and  the  current  step. 

Called  by  MNDO  only. 

PERM  Utility:  Permutes  nl  electrons  of  alpha  or  beta  spin  among 
n2  M.O.s. 

PULAY  Utility:  A  new  converger.  Uses  a  powerful 

mathematical  non-lteratlve  method  for  obtaining  the  SCF  Fock 
matrix.  Principle  Is  that  at  SCF  the  eigenvectors  of  the  Fock 
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and  density  matrices  are  Identical,  so  [F.P]  is  a  measure  of 
the  non-self  consistency.  While  very  powerful,  PULAY  is  not 
universally  applicable.  Used  by  ITER. 

REACT1  Main  sequence:  Uses  reactants  and  products  to  find  the 

transition  state.  A  hypersphere  of  N  dimensions  is  centered  on 
each  moiety,  and  the  radius  steadily  reduced.  The  entity  of 
lower  energy  is  moved,  and  when  the  radius  vanishes,  the 
transition  state  Is  reached.  Called  by  MNDO  only. 

READ  Main  sequence:  Almost  all  the  data  are  read  in  through  READ. 

There  is  a  lot  of  data-checking  in  READ,  but  very  little 
calculation.  Called  by  MNDO. 

READA  Utility:  General  purpose  character  number  reader.  Used  to  enter 
numerical  data  in  the  control  line  as  "  <var1able>-n.nnn  "  where 
<variable>  is  a  mnemonic  such  as  SCFCRT  or  CHARGE. 

Called  by  READ,  FLEPO,  ITER,  FORCE,  and  many  other  subroutines. 

REFER  Utility:  Prints  the  original  references  for  atomic  data. 

If  an  atom  does  not  have  a  reference,  i.e.  it  has  not  been 
parametrized,  then  a  warning  message  will  be  printed  and 
the  calculation  stopped. 

REPP  Utility:  Calculates  the  22  two-electron  reduced  repulsion 
Integrals,  and  the  8  electron-nuclear  attraction  Integrals. 

These  are  in  a  local  coordinate  system.  Arguments:  atomic 
numbers  of  the  two  atoms,  Interatomic  distance,  and  arrays  to 
hold  the  calculated  Integrals.  Called  by  ROTATE  only. 

ROTATE  Utility:  All  the  two-electron  repulsion  Integrals,  the  electron- 
nuclear  attraction  Integrals,  and  the  nuclear-nuclear  repulsion 
term  between  two  atoms  are  calculated  here.  Typically  100  two- 
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electron  Integrals  are  evaluated. 

RSP  Utility:  Rapid  diagonal Izatlon  routine.  Accepts  a  secular 
determinant,  and  produces  a  set  of  eigenvectors  and 
eigenvalues.  The  secular  determinant  Is  destroyed. 

SECOND  Utility:  Contains  IBM-PC  specific  code.  Function  SECOND 
returns  the  number  of  CPU  seconds  elapsed  since  an  arbitrary 
starting  time. 

SPCG  Utility:  Calculates  two-electron  Integral  between  any  four  M.O.'s 
at  the  MNDO  or  MINDO/3  level.  Called  by  MECI  and  WRITE. 

SWAP  Utility:  Used  with  FILL*,  SWAP  ensures  that  a  specified 
M.O.  Is  filled.  Called  by  ITER  only. 

SYMTRY  Utility:  Calculates  values  for  geometric  parameters  from  known 
geometric  parameters  and  symmetry  data.  Called  whenever  GMETRY 
Is  called. 

THERMO  Main  sequence:  After  the  vibrational  frequencies  have  been 
calculated,  THERMO  calculates  thermodynamic  quantities  such  as 
Internal  energy,  heat  capacity,  entropy,  etc,  for  translational, 
vibrational,  and  rotational,  degrees  of  freedom. 

VECPRT  Utility:  Prints  out  a  packed,  lower-half  triangular  matrix. 

The  labeling  of  the  sides  of  the  matrix  depend  on  the  matrix's 
size:  If  It  Is  equal  to  the  number  of  orbitals,  atoms,  or  other. 
Arguments:  The  matrix  to  be  printed,  size  of  matrix.  No  data 
are  changed  by  a  call  of  VECPRT. 

WRITE  Main  sequence:  Most  of  the  results  are  printed  here.  All 

relevant  arrays  are  assumed  to  be  filled.  A  call  of  WRITE  only 
changes  the  number  of  SCF  calls  made,  this  Is  reset  to  zero. 

No  other  data  are  changed.  Called  by  MAIN,  FLEPO,  FORCE. 

XYZINT  Utility:  Converts  from  cartesian  coordinates  Into  Internal. 
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In  order  to  verify  that  MOPAC  Is  working  correctly,  a  large  number 
of  tests  need  to  be  done.  These  take  about  45  minutes  on  a  VAX  11-780, 
and  even  then  the  many  potential  bugs  remain  undetected.  It  Is 
obviously  Impractical  to  ask  users  to  test  MOPAC.  However,  users  must 
be  able  to  verify  the  basic  working  of  MOPAC,  and  to  do  this  the 
following  tests  for  the  elements  have  been  provided. 

Each  element  can  be  tested  by  making  up  a  data-flle  using  estimated 
geometries  and  running  that  file  using  MOPAC.  The  optimized  geometries 
should  give  rise  to  heats  of  formation  as  shown.  Any  difference  greater 
than  0.1  Kcal/mole  Indicates  a  serious  error  In  the  program. 

Caveats 


1.  Geometry  definitions  must  be  correct. 

2.  Heats  of  formation  may  be  too  high  for  certain 
compounds.  This  Is  due  to  a  poor  starting  geometry 
trapping  the  system  In  an  excited  state.  (Affects 
IC1  at  times) 
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ement 

Test  Compound 

Heat  of 
MNDO 

Hydrogen 

CH4 

-11.9 

Lithium 

LI  H 

+23.2 

Beryl  Hum 

BeO 

+38.6 

Boron 

BF3 

-261 .0 

Carbon 

CH4 

-11.9 

N1 trogen 

NH3 

-6.0 

Oxygen 

C02 

-75.0 

Fluorine 

CF4 

-214.2 

Aluminium 

A1F 

-83.6 

Si  1 1  con 

SI  H 

+90.2 

Phosphorus 

PF3 

-229.3 

Sulfur 

H2S 

+  1.7 

Chlorine 

HC1 

-15.3 

Bromine 

HBr 

+3.6 

Formation 

AMI 

-8.8 


-8.8 

-7.3 

-79.8 

-225.7 


-24.6 

-10.5 
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